#Script calculates all statistics presented in Khan et al. (2011) Negative epistasis between beneficial mutations in an evolving bacterial population. This script also contains code to generate all figures in the manuscript and SOM except figures 1 and 4. R version 2.11.0 used in Khan et al. #Formula for calculating epistasis described in da Silva et al. (2010) Fitness epistasis and constraints on adaptation in a human immunodeficiency virus type I protein region. Genetics 184: 293-303. ### Read data folder<-c("~/R/LTmutation_epistasis/Khan_data/") #Folder where input file is and where output files are to be saved data.file<-c("genotype_fitness.txt") #input text file - no spaces or missing values - tab delimited - 1 header row open.file<-paste(folder,data.file,sep="") data=read.table(open.file, header=T) ##Used for other anlyses - retained to keep data appropriately spaced data$genotype=as.character(data$genotype) genes=c("r", "t", "s", "g", "p") ### Functions for reconstructing the mutational network hamming <- function(v1,v2) { u1=strsplit(v1,"")[[1]] u2=strsplit(v2,"")[[1]] return(length(which(u1!=u2))) } benefCount <- function(genotype) { v=as.numeric(strsplit(genotype,"")[[1]]) return(sum(v)) } ### Create binary genotype bin=NULL for (i in 1:length(data$w)) { tmp=strsplit(data$genotype[i],"")[[1]] print(i) print(tmp) bin[i]="" for (j in genes) { bin[i]=paste(bin[i],as.character(as.numeric(is.element(j,tmp))),sep="") } } data$bin=bin ################################################################################## #plot mutation number by fitness - you can uncomment this to plot trajectories that are available - mutation number vs. fitness { fig<-c("figure_S1.pdf") fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=5,height=5) toplot<-aggregate(data$w,by=list(gen=data$genotype,mut=data$mutations),mean) plot(toplot$mut,toplot$x,xlab="mutation number",ylab="relative fitness") #make binary genotypes for the averaged table bina=NULL for (i in 1:length(toplot$x)) { tmp=strsplit(toplot$gen[i],"")[[1]] print(i) print(tmp) bina[i]="" for (j in genes) { bina[i]=paste(bina[i],as.character(as.numeric(is.element(j,tmp))),sep="") } } toplot$bina=bina #transitions with hamming distance==1. for (j in 1:length(toplot$x)) { for (i in (j+1):length(toplot$x)){ genotype1 = toplot$bina[j] genotype2 = toplot$bina[i] if (hamming(genotype1,genotype2)==1) { lines(c(toplot$mut[i],toplot$mut[j]),c(toplot$x[i],toplot$x[j])) } } } toplot=toplot[order(toplot$gen), ] #order to make robust to input order lines(c(0,1,2,3,4,5),c(toplot$x[1],toplot$x[5],toplot$x[13],toplot$x[17],toplot$x[18],toplot$x[19]),lty=1,col="red",lwd=5) #actual trajectory followed in LT evolution experiment dev.off()} ### Sort data according to binary genotype data=data[order(data$bin), ] # Create one column for each locus - identifies rows with corresponding mutation data$locus1=factor(substr(data$bin,1,1)) data$locus2=factor(substr(data$bin,2,2)) data$locus3=factor(substr(data$bin,3,3)) data$locus4=factor(substr(data$bin,4,4)) data$locus5=factor(substr(data$bin,5,5)) geno<-c("anc","r","t","s","g","p","rt","rs","rg","rp","ts","tg","tp","sg","sp","gp","rts","rtg","rtp","rsg","rsp","rgp","tsg","tsp","tgp","sgp","rtsg","rtsp","rtgp","rsgp","tsgp","rtsgp")#list of genotype codes #eqn #1 of da Silva et al. (2010) epidev=numeric(length(geno)) for(i in 1:6){ epidev[i]<-mean(data[data$genotype==geno[i],4]) - mean(data[data$genotype==geno[i],4]) } epidev[7]<-mean(data[data$genotype=="rt",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4]))#rt epidev[8]<-mean(data[data$genotype=="rs",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4]))#rs epidev[9]<-mean(data[data$genotype=="rg",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="g",4]))#rg epidev[10]<-mean(data[data$genotype=="rp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="p",4]))#rp epidev[11]<-mean(data[data$genotype=="ts",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4]))#ts epidev[12]<-mean(data[data$genotype=="tg",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4]))#tg epidev[13]<-mean(data[data$genotype=="tp",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="p",4]))#tp epidev[14]<-mean(data[data$genotype=="sg",4]) - (mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4]))#sg epidev[15]<-mean(data[data$genotype=="sp",4]) - (mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4]))#sp epidev[16]<-mean(data[data$genotype=="gp",4]) - (mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#gp epidev[17]<-mean(data[data$genotype=="rts",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4]))#rts epidev[18]<-mean(data[data$genotype=="rtg",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4]))#rtg epidev[19]<-mean(data[data$genotype=="rtp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="p",4]))#rtp epidev[20]<-mean(data[data$genotype=="rsg",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4]))#rsg epidev[21]<-mean(data[data$genotype=="rsp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4]))#rsp epidev[22]<-mean(data[data$genotype=="rgp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#rgp epidev[23]<-mean(data[data$genotype=="tsg",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4]))#tsg epidev[24]<-mean(data[data$genotype=="tsp",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4]))#tsp epidev[25]<-mean(data[data$genotype=="tgp",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#tgp epidev[26]<-mean(data[data$genotype=="sgp",4]) - (mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#sgp epidev[27]<-mean(data[data$genotype=="rtsg",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4]))#rtsg epidev[28]<-mean(data[data$genotype=="rtsp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4]))#rtsp epidev[29]<-mean(data[data$genotype=="rtgp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#rtgp epidev[30]<-mean(data[data$genotype=="rsgp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#rsgp epidev[31]<-mean(data[data$genotype=="tsgp",4]) - (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#tsgp epidev[32]<-mean(data[data$genotype=="rtsgp",4]) - (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4]))#rtsgp #variance in epistatic deviation equation #2 of da Silva et al. (2010) - this equation is from Goodman 1962 and is a generalization to more than two variables of a formula given in Goodman 1960 and Bohnrstedt 1969 - presented in appendix 1 of 'Genetics and analysis of quantitative traits' M Lynch & B Walsh. Sinauer 1998. temp=numeric(length(geno)) #only the single mutation terms are calculated here - they are put together in the next section! temp2=numeric(length(geno)) for(i in 1:length(geno)){ temp[i]<-(var(data[data$genotype==geno[i],4]) + mean(data[data$genotype==geno[i],4])^2) #first term in equation 2 temp2[i]<-mean(data[data$genotype==geno[i],4])^2 #second term in equation 2 } #variance in epistatic deviation - equation #3 of da Silva et al. (2010) and equation 15.9a (and proof A.13) of Sokal and Rohlf (1995). varepi=numeric(length(geno)) for(i in 1:6){ varepi[i]<-0 } varepi[7]<-var(data[data$genotype=="rt",4]) + (temp[2]*temp[3] - temp2[2]*temp2[3]) #genotypes with >2 mutations - variance is added as stored in temp and temp2 (mut.order as in 'geno' vector) varepi[8]<-var(data[data$genotype=="rs",4]) + (temp[2]*temp[4] - temp2[2]*temp2[4]) #rs varepi[9]<-var(data[data$genotype=="rg",4]) + (temp[2]*temp[5] - temp2[2]*temp2[5])#rg varepi[10]<-var(data[data$genotype=="rp",4]) + (temp[2]*temp[6] - temp2[2]*temp2[6])#rp varepi[11]<-var(data[data$genotype=="ts",4]) + (temp[3]*temp[4] - temp2[3]*temp2[4])#ts varepi[12]<-var(data[data$genotype=="tg",4]) + (temp[3]*temp[5] - temp2[3]*temp2[5])#tg varepi[13]<-var(data[data$genotype=="tp",4]) + (temp[3]*temp[6] - temp2[3]*temp2[6])#tp varepi[14]<-var(data[data$genotype=="sg",4]) + (temp[4]*temp[5] - temp2[4]*temp2[5])#sg varepi[15]<-var(data[data$genotype=="sp",4]) + (temp[4]*temp[6] - temp2[4]*temp2[6])#sp varepi[16]<-var(data[data$genotype=="gp",4]) + (temp[5]*temp[6] - temp2[5]*temp2[6])#gp varepi[17]<-var(data[data$genotype=="rts",4]) + (temp[2]*temp[3]*temp[4] - temp2[2]*temp2[3]*temp2[4])#rts varepi[18]<-var(data[data$genotype=="rtg",4]) + (temp[2]*temp[3]*temp[5] - temp2[2]*temp2[3]*temp2[5])#rtg varepi[19]<-var(data[data$genotype=="rtp",4]) + (temp[2]*temp[3]*temp[6] - temp2[2]*temp2[3]*temp2[6])#rtp varepi[20]<-var(data[data$genotype=="rsg",4]) + (temp[2]*temp[4]*temp[5] - temp2[2]*temp2[4]*temp2[5])#rsg varepi[21]<-var(data[data$genotype=="rsp",4]) + (temp[2]*temp[4]*temp[6] - temp2[2]*temp2[4]*temp2[6])#rsp varepi[22]<-var(data[data$genotype=="rgp",4]) + (temp[2]*temp[5]*temp[6] - temp2[2]*temp2[5]*temp2[6])#rgp varepi[23]<-var(data[data$genotype=="tsg",4]) + (temp[3]*temp[4]*temp[5] - temp2[3]*temp2[4]*temp2[5])#tsg varepi[24]<-var(data[data$genotype=="tsp",4]) + (temp[3]*temp[4]*temp[6] - temp2[3]*temp2[4]*temp2[6])#tsp varepi[25]<-var(data[data$genotype=="tgp",4]) + (temp[3]*temp[5]*temp[6] - temp2[3]*temp2[5]*temp2[6])#tgp varepi[26]<-var(data[data$genotype=="sgp",4]) + (temp[4]*temp[5]*temp[6] - temp2[4]*temp2[5]*temp2[6])#sgp varepi[27]<-var(data[data$genotype=="rtsg",4]) + (temp[2]*temp[3]*temp[4]*temp[5] - temp2[2]*temp2[3]*temp2[4]*temp2[5])#rtsg varepi[28]<-var(data[data$genotype=="rtsp",4]) + (temp[2]*temp[3]*temp[4]*temp[6] - temp2[2]*temp2[3]*temp2[4]*temp2[6])#rtsp varepi[29]<-var(data[data$genotype=="rtgp",4]) + (temp[2]*temp[3]*temp[5]*temp[6] - temp2[2]*temp2[3]*temp2[5]*temp2[6])#rtgp varepi[30]<-var(data[data$genotype=="rsgp",4]) + (temp[2]*temp[4]*temp[5]*temp[6] - temp2[2]*temp2[4]*temp2[5]*temp2[6])#rsgp varepi[31]<-var(data[data$genotype=="tsgp",4]) + (temp[3]*temp[4]*temp[5]*temp[6] - temp2[3]*temp2[4]*temp2[5]*temp2[6])#tsgp varepi[32]<-var(data[data$genotype=="rtsgp",4]) + (temp[2]*temp[3]*temp[4]*temp[5]*temp[6] - temp2[2]*temp2[3]*temp2[4]*temp2[5]*temp2[6])#rtsgp #test for significance of epistatic deviations #a T-test to determine the significance of epistatic interactions - mean from epidev and se = sqrt(varepi)/10 (on average) Pvarepi<-numeric(length(geno)) tvarepi<-numeric(length(geno)) len=NULL for(i in 7:32){ tvarepi[i]<-(epidev[i])/(sqrt(varepi[i])/sqrt(length(data[data$genotype==geno[i],4]))) #'n' calculated for each genotype by last term Pvarepi[i]<-2*pt(-abs(tvarepi[i]),length(data[data$genotype==geno[i],4])) } # two-tailed t-test #net epistatic deviation of higher order interactions - equation #4 of da Silva et al. (2010) Nepidev<-numeric(length(geno)) Nepidev[17] = epidev[17] - (epidev[7]+epidev[8]+epidev[11]) Nepidev[18] = epidev[18] - (epidev[7]+epidev[9]+epidev[12]) Nepidev[19] = epidev[19] - (epidev[7]+epidev[10]+epidev[13]) Nepidev[20] = epidev[20] - (epidev[8]+epidev[9]+epidev[14]) Nepidev[21] = epidev[21] - (epidev[8]+epidev[10]+epidev[15]) Nepidev[22] = epidev[22] - (epidev[9]+epidev[10]+epidev[16]) Nepidev[23] = epidev[23] - (epidev[11]+epidev[12]+epidev[14]) Nepidev[24] = epidev[24] - (epidev[11]+epidev[13]+epidev[15]) Nepidev[25] = epidev[25] - (epidev[12]+epidev[13]+epidev[16]) Nepidev[26] = epidev[26] - (epidev[14]+epidev[15]+epidev[16]) Nepidev[27] = epidev[27] - (epidev[7]+epidev[8]+epidev[9]+epidev[11]+epidev[12]+epidev[14]+epidev[17]+epidev[18]+epidev[20]+epidev[23]) Nepidev[28] = epidev[28] - (epidev[7]+epidev[8]+epidev[10]+epidev[11]+epidev[13]+epidev[15]+epidev[17]+epidev[19]+epidev[21]+epidev[24]) Nepidev[29] = epidev[29] - (epidev[7]+epidev[9]+epidev[10]+epidev[12]+epidev[13]+epidev[16]+epidev[18]+epidev[19]+epidev[25]+epidev[22]) Nepidev[30] = epidev[30] - (epidev[8]+epidev[9]+epidev[10]+epidev[14]+epidev[15]+epidev[16]+epidev[20]+epidev[22]+epidev[26]+epidev[21]) Nepidev[31] = epidev[31] - (epidev[11]+epidev[12]+epidev[13]+epidev[14]+epidev[15]+epidev[16]+epidev[23]+epidev[24]+epidev[25]+epidev[26]) Nepidev[32] = epidev[32] - (epidev[7]+epidev[8]+epidev[9]+epidev[10]+epidev[11]+epidev[12]+epidev[13]+epidev[14]+epidev[15]+epidev[16]+epidev[17]+epidev[18]+epidev[19]+epidev[20]+epidev[21]+epidev[22]+epidev[23]+epidev[24]+epidev[25]+epidev[26]+epidev[27]+epidev[28]+epidev[29]+epidev[30]+epidev[31]) #variance in net epistatic deviation - equation #6 of da Silva et al. (2010) varNepidev<-numeric(length(geno)) varNepidev[17] = varepi[17] + varepi[7] + varepi[8] + varepi[11] varNepidev[18] = varepi[18] + varepi[7]+ varepi[9]+ varepi[12] varNepidev[19] = varepi[19] + varepi[7]+ varepi[10]+ varepi[13] varNepidev[20] = varepi[20] + varepi[8]+ varepi[9]+ varepi[14] varNepidev[21] = varepi[21] + varepi[8]+ varepi[10]+ varepi[15] varNepidev[22] = varepi[22] + varepi[9]+ varepi[10]+ varepi[16] varNepidev[23] = varepi[23] + varepi[11]+ varepi[12]+ varepi[14] varNepidev[24] = varepi[24] + varepi[11]+ varepi[13]+ varepi[15] varNepidev[25] = varepi[25] + varepi[12]+ varepi[13]+ varepi[16] varNepidev[26] = varepi[26] + varepi[14]+ varepi[15]+ varepi[16] varNepidev[27] = varepi[27] + varepi[7]+ varepi[8]+ varepi[9]+ varepi[11]+ varepi[12]+ varepi[14]+ varepi[17]+ varepi[18]+ varepi[20]+ varepi[23] varNepidev[28] = varepi[28] + varepi[7]+ varepi[8]+ varepi[10]+ varepi[11]+ varepi[13]+ varepi[15]+ varepi[17]+ varepi[19]+ varepi[21]+ varepi[24] varNepidev[29] = varepi[29] + varepi[7]+ varepi[9]+ varepi[10]+ varepi[12]+ varepi[13]+ varepi[16]+ varepi[18]+ varepi[19]+ varepi[25]+ varepi[22] varNepidev[30] = varepi[30] + varepi[8]+ varepi[9]+ varepi[10]+ varepi[14]+ varepi[15]+ varepi[16]+ varepi[20]+ varepi[22]+ varepi[26]+ varepi[21] varNepidev[31] = varepi[31] + varepi[11]+ varepi[12]+ varepi[13]+ varepi[14]+ varepi[15]+ varepi[16]+ varepi[23]+ varepi[24]+ varepi[25]+ varepi[26] varNepidev[32] = varepi[32] + varepi[7]+varepi[8]+varepi[9]+varepi[10]+varepi[11]+varepi[12]+varepi[13]+varepi[14]+varepi[15]+varepi[16]+varepi[17]+varepi[18]+varepi[19]+varepi[20]+varepi[21]+varepi[22]+varepi[23]+varepi[24]+varepi[25]+varepi[26]+varepi[27]+varepi[28]+varepi[29]+varepi[30]+varepi[31] #a T-test to determine the significance of epistatic interactions - mean from Nepidev and se = sqrt(varNepidev)/22 (on average) PNepidev<-numeric(length(geno)) tNepidev<-numeric(length(geno)) for(i in 17:32){ tNepidev[i]<-(Nepidev[i])/(sqrt(varNepidev[i])/sqrt(length(data[data$genotype==geno[i],4]))) PNepidev[i]<-2*pt(-abs(tNepidev[i]),length(data[data$genotype==geno[i],4]))} #magnitude of epistasis equation #8 of da Silva et al. (2010) - gives overall (not net!) effect of epistasis for each genotype magepi<-numeric(length(geno)) magepi[7]= log(mean(data[data$genotype=="rt",4])/(mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])),10)#rt magepi[8]<-log(mean(data[data$genotype=="rs",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])),10)#rs magepi[9]<-log(mean(data[data$genotype=="rg",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="g",4])),10)#rg magepi[10]<-log(mean(data[data$genotype=="rp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="p",4])),10)#rp magepi[11]<-log(mean(data[data$genotype=="ts",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])),10)#ts magepi[12]<-log(mean(data[data$genotype=="tg",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4])),10)#tg magepi[13]<-log(mean(data[data$genotype=="tp",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="p",4])),10)#tp magepi[14]<-log(mean(data[data$genotype=="sg",4]) / (mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])),10)#sg magepi[15]<-log(mean(data[data$genotype=="sp",4]) / (mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4])),10)#sp magepi[16]<-log(mean(data[data$genotype=="gp",4]) / (mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#gp magepi[17]<-log(mean(data[data$genotype=="rts",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])),10)#rts magepi[18]<-log(mean(data[data$genotype=="rtg",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4])),10)#rtg magepi[19]<-log(mean(data[data$genotype=="rtp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="p",4])),10)#rtp magepi[20]<-log(mean(data[data$genotype=="rsg",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])),10)#rsg magepi[21]<-log(mean(data[data$genotype=="rsp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4])),10)#rsp magepi[22]<-log(mean(data[data$genotype=="rgp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#rgp magepi[23]<-log(mean(data[data$genotype=="tsg",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])),10)#tsg magepi[24]<-log(mean(data[data$genotype=="tsp",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4])),10)#tsp magepi[25]<-log(mean(data[data$genotype=="tgp",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#tgp magepi[26]<-log(mean(data[data$genotype=="sgp",4]) / (mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#sgp magepi[27]<-log(mean(data[data$genotype=="rtsg",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])),10)#rtsg magepi[28]<-log(mean(data[data$genotype=="rtsp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="p",4])),10)#rtsp magepi[29]<-log(mean(data[data$genotype=="rtgp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#rtgp magepi[30]<-log(mean(data[data$genotype=="rsgp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#rsgp magepi[31]<-log(mean(data[data$genotype=="tsgp",4]) / (mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#tsgp magepi[32]<-log(mean(data[data$genotype=="rtsgp",4]) / (mean(data[data$genotype=="r",4])*mean(data[data$genotype=="t",4])*mean(data[data$genotype=="s",4])*mean(data[data$genotype=="g",4])*mean(data[data$genotype=="p",4])),10)#rtsgp final<-NULL finalNum<-NULL fit<-c(mean(data[data$genotype=="anc",4]),mean(data[data$genotype=="r",4]),mean(data[data$genotype=="t",4]),mean(data[data$genotype=="s",4]),mean(data[data$genotype=="g",4]),mean(data[data$genotype=="p",4]),mean(data[data$genotype=="rt",4]),mean(data[data$genotype=="rs",4]),mean(data[data$genotype=="rg",4]),mean(data[data$genotype=="rp",4]),mean(data[data$genotype=="ts",4]),mean(data[data$genotype=="tg",4]),mean(data[data$genotype=="tp",4]),mean(data[data$genotype=="sg",4]),mean(data[data$genotype=="sp",4]),mean(data[data$genotype=="gp",4]),mean(data[data$genotype=="rts",4]),mean(data[data$genotype=="rtg",4]),mean(data[data$genotype=="rtp",4]),mean(data[data$genotype=="rsg",4]),mean(data[data$genotype=="rsp",4]),mean(data[data$genotype=="rgp",4]),mean(data[data$genotype=="tsg",4]),mean(data[data$genotype=="tsp",4]),mean(data[data$genotype=="tgp",4]),mean(data[data$genotype=="sgp",4]),mean(data[data$genotype=="rtsg",4]),mean(data[data$genotype=="rtsp",4]),mean(data[data$genotype=="rtgp",4]),mean(data[data$genotype=="rsgp",4]),mean(data[data$genotype=="tsgp",4]),mean(data[data$genotype=="rtsgp",4])) #Combine all calculations from above to create a summary table final3<-data.frame(geno,fit,epidev,varepi,tvarepi,Pvarepi,Nepidev,varNepidev,tNepidev,PNepidev,magepi) #predicted (multiplicative) fitness - i.e. null in absence of epistasis pred_fit=c(fit[1],fit[2],fit[3],fit[4],fit[5],fit[6],fit[2]*fit[3],fit[2]*fit[4],fit[2]*fit[5],fit[2]*fit[6], fit[3]*fit[4],fit[3]*fit[5],fit[3]*fit[6],fit[4]*fit[5],fit[4]*fit[6],fit[5]*fit[6], fit[2]*fit[3]*fit[4],fit[2]*fit[3]*fit[5],fit[2]*fit[3]*fit[6],fit[2]*fit[4]*fit[5],fit[2]*fit[4]*fit[6], fit[2]*fit[5]*fit[6],fit[3]*fit[4]*fit[5],fit[3]*fit[4]*fit[6],fit[3]*fit[5]*fit[6],fit[4]*fit[5]*fit[6],fit[2]*fit[3]*fit[4]*fit[5],fit[2]*fit[3]*fit[4]*fit[6],fit[2]*fit[3]*fit[5]*fit[6],fit[2]*fit[4]*fit[5]*fit[6],fit[3]*fit[4]*fit[5]*fit[6],fit[2]*fit[3]*fit[4]*fit[5]*fit[6]) #key for mutation presence/absence Fgene=final3[7:32,] rbs<- c(1,1,1,1,0,0,0,0,0,0,1,1,1,1,1,1,0,0,0,0,1,1,1,1,0,1) topA<- c(1,0,0,0,1,1,1,0,0,0,1,1,1,0,0,0,1,1,1,0,1,1,1,0,1,1) spoT<- c(0,1,0,0,1,0,0,1,1,0,1,0,0,1,1,0,1,1,0,1,1,1,0,1,1,1) glmUS<-c(0,0,1,0,0,1,0,1,0,1,0,1,0,1,0,1,1,0,1,1,1,0,1,1,1,1) pykF<- c(0,0,0,1,0,0,1,0,1,1,0,0,1,0,1,1,0,1,1,1,0,1,1,1,1,1) Fgene <-data.frame(Fgene,rbs,topA,spoT,glmUS,pykF,pred_fit[7:32]) attach(Fgene) #95% CI of mean genotype fitness values ci=NULL for(i in 1:32){ meanGen=mean(data[data$genotype==geno[i],4]) SE<-sd(data[data$genotype==geno[i],4])/sqrt(length(data[data$genotype==geno[i],4])) tStat<-abs(meanGen-1)/SE temp<-SE*qt(.975,df=34) ci<-c(ci,temp) } ###Table S2###### tableS2<-data.frame(final3[,2],ci,final3[,3],final3[,6],final3[,11]) tableS2<-round(tableS2,3) #rounded values as presented in SOM tableS2<-cbind(final3[,1],tableS2) names(tableS2)=c("genotype","w","95ci","Abs_epistasis","P","Rel_epistasis") ################histograms of epistasis effects############################ { fig<-c("Fig_2.pdf") fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=4,height=3.5) par(fig=c(0,0.55,0,1)) par(mar=c(5,3.8,2,0)) hist(epidev[7:32],xlab="absolute epistasis",ylab="count",main="",ylim=c(0,6),col="gray",breaks=10,freq=TRUE,cex.axis=0.7,cex.lab=1,xaxp=c(-0.15,0.15,2),yaxp=c(0,6,2),xlim=c(-0.15,0.15)) #7:32 because the single mutants don't have epistasis calculated lines(c(mean(epidev[7:32]),mean(epidev[7:32])),c(-1,6),col="RED") text(-0.1,5.9,"A",cex=1,font=2) par(fig=c(0.55,1,0,1),new=T) par(mar=c(5,0.7,2,1)) hist(magepi[7:32],xlab="relative epistasis",ylab="",main="",ylim=c(0,6),col="gray",breaks=15,freq=TRUE,axes=FALSE,cex.lab=1,xaxp=c(-0.05,0.05,2),xlim=c(-0.05,0.05)) #7:32 because the single mutants don't have epistasis calculated #axis(2,labels=FALSE) #x-axis axis(1,at=c(-0.05,0,0.05),cex.axis=.7) lines(c(mean(magepi[7:32]),mean(magepi[17:32])),c(-1,6),col="RED") text(-0.03,5.9,"B",cex=1,font=2) dev.off() } #####################separating genotypes with and without pykF ################################ fig<-c("figure_s4.pdf") fig.file<-paste(folder,fig,sep="") superhist2pdf <- function(x, filename = fig.file, dev = "pdf", title = "", nbreaks ="Sturges") { junk = NULL grouping = NULL for(i in 1:length(x)) { junk = c(junk,x[[i]]) grouping <- c(grouping, rep(i,length(x[[i]]))) } grouping <- factor(grouping) n.gr <- length(table(grouping)) xr <- range(junk) histL <- tapply(junk, grouping, hist, breaks=nbreaks, plot = FALSE) maxC <- max(sapply(lapply(histL, "[[", "counts"), max)) if(dev == "pdf") { pdf(filename, version = "1.4") } else{} if((TC <- transparent.cols <- .Device %in% c("pdf", "png"))) { cols <- hcl(h = seq(240, by=120 , length = n.gr), l = 60, alpha = 0.5) }#240+120 give initial and offset of color in the HCL 360 description else { h.den <- c(10, 15, 20) h.ang <- c(45, 15, -30) } if(TC) { plot(histL[[1]], xlim = c(-0.06,0.06), ylim= c(0, 6), col = cols[1], xlab = "relative epistasis", ylab="count",main = title) } if(!transparent.cols) { for(j in 2:n.gr) plot(histL[[j]], add = TRUE, density = h.den[j], angle = h.ang[j]) } else { for(j in 2:n.gr) plot(histL[[j]], add = TRUE, col = cols[j]) } invisible() ###mean and 95% CI of genotype w and w/o pykF test<-t.test((magepi[7:32][Fgene[,16]=="0"])) testp<-t.test((magepi[7:32][Fgene[,16]=="1"])) test2<-test[[4]] test2p<-testp[[4]] lines(c(test2[1],test2[2]),c(5.5,5.5),col="Blue") lines(c(test2p[1],test2p[2]),c(5.5,5.5),col="red") points(mean(magepi[7:32][Fgene[,16]=="0"]),5.5,col="blue",pch=19,cex=1.2) points(mean(magepi[7:32][Fgene[,16]=="1"]),5.5,col="red",pch=19,cex=1.2) if( dev == "pdf") { dev.off() } } # Call histogram function d1 = magepi[7:32][Fgene[,16]=="0"] d2 = magepi[7:32][Fgene[,16]=="1"] l1 = list(d1,d2) # the input object MUST be a list! superhist2pdf(l1, nbreaks="Sturges") ####Alternative figure_S4 - side by side and stacked histograms d1 = magepi[7:32][Fgene[,16]=="0"] d2 = magepi[7:32][Fgene[,16]=="1"] h<-hist(d1,breaks=c(-.05,-.04,-.03,-.02,-.01,0,.01,.02,.03,.04,.05),plot=FALSE) #creates counts in each bin to be plotted using 'barplot' h2<-hist(d2,breaks=c(-.05,-.04,-.03,-.02,-.01,0,.01,.02,.03,.04,.05),plot=FALSE) h3<-rbind(h$counts,h2$counts) fig<-c("figure_s4_alt.pdf") fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=4,height=4) plot(1,1,col="white",axes=FALSE,ylab="count",xlim=c(-0.05,0.05), ylim=c(0,6),xlab="relative epistasis") #empty plot to set scale for axes axis(1,at=c(-0.05,0,0.05)) test<-t.test((magepi[7:32][Fgene[,16]=="0"])) testp<-t.test((magepi[7:32][Fgene[,16]=="1"])) test2<-test[[4]] test2p<-testp[[4]] lines(c(test2[1],test2[2]),c(6,6),col="Blue") lines(c(test2p[1],test2p[2]),c(6,6),col="red") points(mean(magepi[7:32][Fgene[,16]=="0"]),6,col="blue",pch=19,cex=1.2) points(mean(magepi[7:32][Fgene[,16]=="1"]),6,col="red",pch=19,cex=1.2) par(fig=c(0,1,0,1),new=T) par(mar=c(5.5,4,3.5,2)) barplot(as.matrix(h3),horiz=F,col=c("blue","red"),beside=TRUE,ylim=c(0,6)) #side by side #barplot(as.matrix(h3),space=0,horiz=F,col=c("blue","red"),ylim=c(0,6)) #stacked dev.off() ############################################################################################ #################-----------------------Figure Function-------------######################## fig.func<-function(dep,indep,mut,xla,yla,mutBin){ #ANCOVA and correlations ancova<-lm(dep[7:32]~indep[7:32]*mutBin) #full model (only one mutation) column11 = magepi ancova1<-update(ancova, ~ . -indep[7:32]:mutBin) #remove interaction term - leave 1 slope and 2 intercepts ancova2<-update(ancova1, ~ . -mutBin) #remove mutation - leaving only initial fitness - 1 slope and 1 intercept model<-anova(ancova,ancova1) model1<-anova(ancova1,ancova2) corr<-cor.test(dep[7:32],indep[7:32]) corrM<-cor.test(dep[7:32][mutBin=="1"],indep[7:32][mutBin=="1"]) corrM2<-cor.test(dep[7:32][mutBin=="0"],indep[7:32][mutBin=="0"]) #plot plot(indep[7:32],dep[7:32],pch=16+as.numeric(mutBin),col=ifelse(mutBin=="0","blue","red"),xlab=xla,ylab=yla, xlim=c(1,1.35),ylim=c(-0.04,0.04)) #plots shape and color depending on focal mutation abline(0,0,lty=2,col="gray") abline(lm(dep[7:32][mutBin=="0"]~indep[7:32][mutBin=="0"]),lty=2,col="blue") abline(lm(dep[7:32][mutBin=="1"]~indep[7:32][mutBin=="1"]),lty=2,col="red") abline(lm(dep[7:32]~indep[7:32]),lty=2) #correlation P values - all values and separately with and without focal mutation text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.985),mut,col="red") text(min(indep)+((max(indep)-min(indep))*.83),min(dep)+((max(dep)-min(dep))*.825),"r = ",cex=.8,col="red") ##If statements are if p values are plotted - not necessary if plotting r # if(corrM[[4]] > 0.001) {text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.825),round(corrM[[4]],digits=3),cex=.8,col="red")} # else{text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.825),"< 0.001",cex=.8,col="red")} text(min(indep)+((max(indep)-min(indep))*.83),min(dep)+((max(dep)-min(dep))*.875),"r = ",cex=.8,col="blue") # if(corrM2[[3]] > 0.001) {text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.875),round(corrM2[[4]],digits=4),cex=.8,col="blue")} # else{text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.875),"< 0.001",cex=.8,col="blue")} text(min(indep)+((max(indep)-min(indep))*.83),min(dep)+((max(dep)-min(dep))*.925),"r = ",cex=.8) # if(corr[[3]] > 0.001) {text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.925),round(corr[[4]],digits=4),cex=.8,)} # else{text(min(indep)+((max(indep)-min(indep))*.9),min(dep)+((max(dep)-min(dep))*.925),"< 0.001",cex=.8,)} #mutation and interaction term P values text(min(indep)+((max(indep)-min(indep))*.3),min(dep)+((max(dep)-min(dep))*.12),"-interaction",cex=.8) if(model[[2,6]] > 0.001){text(min(indep)+((max(indep)-min(indep))*.3),min(dep)+((max(dep)-min(dep))*.05),round(model[[2,6]],digits=3),cex=.8)} else{text(min(indep)+((max(indep)-min(indep))*.3),min(dep)+((max(dep)-min(dep))*.05),"<0.001",cex=.8)} text(min(indep)+((max(indep)-min(indep))*.8),min(dep)+((max(dep)-min(dep))*.12),"-mutation",cex=.8) if(model1[[2,6]] > 0.001){text(min(indep)+((max(indep)-min(indep))*.8),min(dep)+((max(dep)-min(dep))*.05),round(model1[[2,6]],digits=3),cex=.8)} else{text(min(indep)+((max(indep)-min(indep))*.8),min(dep)+((max(dep)-min(dep))*.05),"<0.001",cex=.8)} } ############################################################################################## ############################################################################################### ###Each section calls fig.func to generate set of 5 fig - 1 for each focal mutation fig<-c("figure_S3.pdf") fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=8.5,height=11) par(mfrow = c(3,2)) fig.func(magepi,pred_fit,"+rbs-ev","expected multiplicative fitness","relative epistasis",Fgene[,12]) #Fgene gives a binary code for presence/absence of focal mutation fig.func(magepi,pred_fit,"+topA-ev","expected multiplicative fitness","relative epistasis",Fgene[,13]) fig.func(magepi,pred_fit,"+spoT-ev","expected multiplicative fitness","relative epistasis",Fgene[,14]) fig.func(magepi,pred_fit,"+glmUS-ev","expected multiplicative fitness","relative epistasis",Fgene[,15]) fig.func(magepi,pred_fit,"+pykF-ev","expected multiplicative fitness","relative epistasis",Fgene[,16]) dev.off() {fig<-c("Fig_3.pdf") fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=5,height=5) plot(pred_fit[7:32],magepi[7:32],pch=16+as.numeric(Fgene[,16]),col=ifelse(Fgene[,16]=="0","blue","red"),xlim=c(0.95,1.35),yaxp=c(-0.04,0.04,2),ylim=c(-0.04,0.04),xlab="expected multiplicative fitness",ylab="relative epistasis",cex.lab=1.4) #plots shape and color depending on presence of focal mutation abline(0,0,lty=2,col="gray") abline(lm(magepi[7:32][Fgene[,16]=="0"]~pred_fit[7:32][Fgene[,16]=="0"]),lty=2,col="blue") abline(lm(magepi[7:32][Fgene[,16]=="1"]~pred_fit[7:32][Fgene[,16]=="1"]),lty=2,col="red") dev.off()} {fig<-c("figure_S5.pdf") #Absolute epistasis vs. expected fitness fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=5,height=5) plot(pred_fit[7:32],epidev[7:32],pch=16+as.numeric(Fgene[,16]),col=ifelse(Fgene[,16]=="0","blue","red"),xlim=c(0.95,1.35),ylim=c(-0.15,0.15),yaxp=c(-.15,.15,2),xlab="expected multiplicative fitness",ylab="absolute epistasis") #plots shape and color depending on presence of focal mutation abline(0,0,lty=2,col="gray") abline(lm(epidev[7:32][Fgene[,16]=="0"]~pred_fit[7:32][Fgene[,16]=="0"]),lty=2,col="blue") abline(lm(epidev[7:32][Fgene[,16]=="1"]~pred_fit[7:32][Fgene[,16]=="1"]),lty=2,col="red") dev.off()} {fig<-c("figure_S2.pdf") #Absolute epistasis vs. expected fitness fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=5,height=5) plot(pred_fit[7:32],fit[7:32],pch=16,xlim=c(0.95,1.35),ylim=c(0.95,1.35),xlab="expected relative fitness",ylab="observed relative fitness") #plots shape and color depending on presence of focal mutation abline(0,1,lty=2,col="gray") dev.off()} {fig<-c("figure_S6.pdf") #Absolute epistasis vs. expected fitness fig.file<-paste(folder,fig,sep="") pdf(fig.file,onefile=F,width=5,height=5) plot(pred_fit[7:32],fit[7:32],pch=16+as.numeric(Fgene[,16]),col=ifelse(Fgene[,16]=="0","blue","red"),xlim=c(0.95,1.35),ylim=c(0.95,1.35),xlab="expected multiplicative fitness",ylab="observed fitness") #plots shape and color depending on presence of focal mutation abline(0,1,lty=2,col="gray") ##These lines calculated from Major axis regression -- below for details abline(0.4,0.628,col="blue") abline(0.26,0.748,col="blue",lty=2) abline(0.532,0.52,,col="blue",lty=2) abline(0.213,0.85,col="red") abline(0.0064,1.028,col="red",lty=2) abline(0.389,0.698,col="red",lty=2) dev.off()} #########Statistics reported in manuscript text #Main text: t.test of absolute and relative epistasis #absolute epistasis t.test(epidev[7:32]) #relative epistasis t.test(magepi[7:32]) #Main text: t.test - one significant decrease in fitness without multiple correction t.test(data[data$genotype=="t",4],data[data$genotype=="rt",4]) #P=009478>P=0.05/80 tests for Bonferroni multiple correction #Main text: correlation between epistasis and expected fitness #absolute epistasis vs. expected fitnes cor.test(pred_fit[7:32],epidev[7:32]) #relative epistasis vs. expected fitness cor.test(pred_fit[7:32],magepi[7:32]) #Main text: t.test of epistasis with and without pykF #relative epistasis of genotypes not having the pykF-ev mutation t.test((magepi[7:32][Fgene[,16]=="0"])) #relative epistasis of genotypes that have the pykF-ev mutation t.test((magepi[7:32][Fgene[,16]=="1"])) #Main text: correlation between relative epistasis and expected fitness -- genotypes with and without pykF #relative epistasis vs. expected fitnes -- genotypes without pykF cor.test(pred_fit[7:32][Fgene[,16]=="0"],magepi[7:32][Fgene[,16]=="0"]) #relative epistasis vs. expected fitnes -- genotypes with pykF cor.test(pred_fit[7:32][Fgene[,16]=="1"],magepi[7:32][Fgene[,16]=="1"]) ##Main text: major axis regression -- slope of expectd vs. observed genotype fitness w/wo pykF mutaiton library(lmodel2) #genotypes without pykF. Results from here used to calculate slopes and 95%CI as plotted in figure S6. lmodel2(fit[7:32][Fgene[,16]=="0"]~pred_fit[7:32][Fgene[,16]=="0"], nperm=99) #genotypes with pykF. Results from here used to calculate slopes and 95%CI as plotted in figure S6. lmodel2(fit[7:32][Fgene[,16]=="1"]~pred_fit[7:32][Fgene[,16]=="1"], nperm=99) #SOM table S2 #Note: The rounding function of R converts small negative absolute and relative epistasis #values for the 'rs' genotype and the relative epistasis of 'rtsp' to '0.000'. In these cases, #the original negative sign is preserved in SOM for consistency with plotted data - e.g. figure S4 tableS2 ################################################################## #table_s4 ANCOVA summarized in figure S3 and all 2-way interactions epi<-magepi[7:32] fit<-pred_fit[7:32] rbs<-Fgene[,12] top<-Fgene[,13] spo<-Fgene[,14] glm<-Fgene[,15] pyk<-Fgene[,16] #rbs focal ancova<-aov(epi~fit*rbs*top) ancova1<-update(ancova, ~ . -fit:rbs:top) ancova2<-update(ancova1, ~ . -rbs:top) anova(ancova,ancova2) ancova<-aov(epi~fit*rbs*spo) ancova1<-update(ancova, ~ . -fit:rbs:spo) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -rbs:spo) anova(ancova,ancova2) ancova<-aov(epi~fit*rbs*glm) ancova1<-update(ancova, ~ . -fit:rbs:glm) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -rbs:glm) anova(ancova,ancova2) ancova3<-aov(epi~fit*rbs*pyk) ancova4<-update(ancova3, ~ . -fit:rbs:pyk) # 2-way interactions including focal mutation ancova5<-update(ancova4, ~ . -rbs:pyk) anova(ancova3,ancova5) ancova<-aov(epi~fit*rbs) ancova1<-update(ancova, ~ . -fit:rbs) #remove interaction term ancova2<-update(ancova1, ~ . -rbs) #remove mutation main effect anova(ancova,ancova1) anova(ancova1,ancova2) #topA focal ancova<-aov(epi~fit*top*rbs) ancova1<-update(ancova, ~ . -fit:top:rbs) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -top:rbs) anova(ancova,ancova2) ancova<-aov(epi~fit*top*spo) ancova1<-update(ancova, ~ . -fit:top:spo) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -top:spo) anova(ancova,ancova2) ancova<-aov(epi~fit*top*glm) ancova1<-update(ancova, ~ . -fit:top:glm) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -top:glm) anova(ancova,ancova2) ancova<-aov(epi~fit*top*pyk) ancova1<-update(ancova, ~ . -fit:top:pyk) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -top:pyk) anova(ancova,ancova2) ancova<-aov(epi~fit*top) ancova1<-update(ancova, ~ . -fit:top) #remove interaction term ancova2<-update(ancova1, ~ . -top) #remove mutation main effect anova(ancova,ancova1) anova(ancova1,ancova2) #spoT focal ancova<-aov(epi~fit*spo*rbs) ancova1<-update(ancova, ~ . -fit:spo:rbs) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -spo:rbs) anova(ancova,ancova2) ancova<-aov(epi~fit*spo*top) ancova1<-update(ancova, ~ . -fit:spo:top) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -spo:top) anova(ancova,ancova2) ancova<-aov(epi~fit*spo*glm) ancova1<-update(ancova, ~ . -fit:spo:glm) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -spo:glm) anova(ancova,ancova2) ancova<-aov(epi~fit*spo*pyk) ancova1<-update(ancova, ~ . -fit:spo:pyk) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -spo:pyk) anova(ancova,ancova2) ancova<-aov(epi~fit*spo) ancova1<-update(ancova, ~ . -fit:spo) #remove interaction term ancova2<-update(ancova1, ~ . -spo) #remove mutation main effect anova(ancova,ancova1) anova(ancova1,ancova2) #glmUS focal ancova<-aov(epi~fit*glm*rbs) ancova1<-update(ancova, ~ . -fit:glm:rbs) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -glm:rbs) anova(ancova,ancova2) ancova<-aov(epi~fit*glm*top) ancova1<-update(ancova, ~ . -fit:glm:top) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -glm:top) anova(ancova,ancova2) ancova<-aov(epi~fit*glm*spo) ancova1<-update(ancova, ~ . -fit:glm:spo) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -glm:spo) anova(ancova,ancova2) ancova<-aov(epi~fit*glm*pyk) ancova1<-update(ancova, ~ . -fit:glm:pyk) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -glm:pyk) anova(ancova,ancova2) ancova<-aov(epi~fit*glm) ancova1<-update(ancova, ~ . -fit:glm) #remove interaction term ancova2<-update(ancova1, ~ . -glm) #remove mutation main effect anova(ancova,ancova1) anova(ancova1,ancova2) #pykF focal ancova<-aov(epi~fit*pyk*rbs) ancova1<-update(ancova, ~ . -fit:pyk:rbs) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -pyk:rbs) anova(ancova,ancova2) ancova<-aov(epi~fit*pyk*top) ancova1<-update(ancova, ~ . -fit:pyk:top) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -pyk:top) anova(ancova,ancova2) ancova<-aov(epi~fit*pyk*spo) ancova1<-update(ancova, ~ . -fit:pyk:spo) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -pyk:spo) anova(ancova,ancova2) ancova<-aov(epi~fit*pyk*glm) ancova1<-update(ancova, ~ . -fit:pyk:glm) # 2-way interactions including focal mutation ancova2<-update(ancova1, ~ . -pyk:glm) anova(ancova,ancova2) ancova<-aov(epi~fit*pyk) ancova1<-update(ancova, ~ . -fit:pyk) #remove interaction term ancova2<-update(ancova1, ~ . -pyk) #remove mutation main effect anova(ancova,ancova1) anova(ancova1,ancova2)