Distinct evolutionary trajectories of two integration centres, the central complex and mushroom bodies, across Heliconiini butterflies.

Code ↔ Paper

3 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 3 matches

[1] § Materials and methods › Image processing, 3D segmentation, and annotation for fine-scale analysis ↔ Heliconiini-CX_SCRIPT.R, lines 2226–2307 · score 0.63 · bulb volume, ER neuron, metrics, quantified, segment, cell
[2] § Results › GABA-ergic ER neurons increase in number, with increased innervation of the fan-shaped body in pollen feeders ↔ Heliconiini-CX_SCRIPT.R, lines 2226–2307 · score 0.63 · predict bulb, bulb volume, ER neuron, segmented, species, sex
[3] § Results › Isolated impact of mushroom body expansion within the central complex, AOTU, and POTU circuit ↔ Heliconiini-CX_SCRIPT.R, lines 1–44 · score 0.54 · closer inspection, nested model, simplest, MCMCglmm, rCBR, sex

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The paper is loaded by your browser, with JavaScript: read it at doi.org or on Europe PMC.

The authors' code

R · 2,308 lines · 87 KB · CC-BY-4.0 · 3 matches

##### PREREQUISITES #####
rm(list=ls())
library(tidyverse)
library(ggpubr)
library(phangorn)
library(MCMCglmm)
library(phytools)
library(bbmle)
library(MuMIn)
library(plotly)
library(ggbeeswarm)
library(wesanderson)
library(ggResidpanel)
library(emmeans)
##### #SEX DIFFERENCES (FIGURE S1) #####
setwd("C:/R-Analysis")
cx.data=read.table(file="Hel-CX-MB_vs4.txt", sep="\t", header=TRUE)
str(cx.data)
#test whether pervasive sex differences exist that need closer inspection.
#perform three different nested models, then compare their fit, to use the best fit and simplest to interpret.
#TREE
tree = read.nexus("Heliconiini.trees")
tree=force.ultrametric(tree,method="nnls")
plot(tree)
inv.phylo=inverseA(tree,nodes="TIPS",scale=TRUE)
##MCMCglmm
# set priors
prior=list(G=list(G1=list(V=1,nu=0.02)),R=list(V=1,nu=0.02)) #standard priors as explored by SHM
## define models
model.sexA=log10(Total.CX.woNO)~log10(rCBR)+SEX.mod
## run model
model.sexA.1=MCMCglmm(model.sexA, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.sexA.2=MCMCglmm(model.sexA, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.sexA.1$Sol,model.sexA.2$Sol))
#good The test statistic is called the scale reduction factor. The closer this factor is to 1, the better the convergence of our chains. In practice, values below 1.1 can be acceptable and values below 1.02 are good.
gelman.diag(mcmc.list(model.sexA.1$VCV,model.sexA.2$VCV))
plot(mcmc.list(model.sexA.1$VCV,model.sexA.2$VCV))
## autocorrelation
autocorr(model.sexA.1$Sol)
autocorr(model.sexA.1$VCV)
#alternative model B
## define models
model.sexB=log10(Total.CX.woNO)~log10(rCBR)+SEX.mod+HelNonHel
## run model
model.sexB.1=MCMCglmm(model.sexB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.sexB.2=MCMCglmm(model.sexB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.sexB.1$Sol,model.sexB.2$Sol))
#okay
gelman.diag(mcmc.list(model.sexB.1$VCV,model.sexB.2$VCV))
#okay? phylo a bit
plot(mcmc.list(model.sexB.1$VCV,model.sexB.2$VCV))
## autocorrelation
autocorr(model.sexB.1$Sol)
autocorr(model.sexB.1$VCV)
#alternative model C
## define models
model.sexC=log10(Total.CX.woNO)~log10(rCBR)+SEX.mod*HelNonHel
## run model
model.sexC.1=MCMCglmm(model.sexC, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.sexC.2=MCMCglmm(model.sexC, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.sexC.1$Sol,model.sexC.2$Sol))
gelman.diag(mcmc.list(model.sexC.1$VCV,model.sexC.2$VCV))
plot(mcmc.list(model.sexC.1$VCV,model.sexC.2$VCV))
## autocorrelation
autocorr(model.sexC.1$Sol)
autocorr(model.sexC.1$VCV)
#comparison of nested models.
DIC(model.sexA.1)
DIC(model.sexB.1)
DIC(model.sexC.1)
#differences are miniscule at best, hence its worth checking the model A still as its the simplest by far and performs similarly.
summary(model.sexA.1)
#p. val is indeed 0.031
#plot
log.rCBR=log10(cx.data$rCBR)
log.CX=log10(cx.data$Total.CX.woNO)
log.CBU=log10(cx.data$CBU)
log.CBL=log10(cx.data$CBL)
log.PB=log10(cx.data$PB)
log.AOTU=log10(cx.data$AOTU)
log.POTU=log10(cx.data$POTU)
#add this to cx.data
cx.data.log=cx.data %>%
mutate(log.rCBR = log.rCBR,log.CX=log.CX,log.CBU=log.CBU,log.CBL=log.CBL,log.PB=log.PB,log.AOTU=log.AOTU,log.POTU=log.POTU)
str(cx.data.log)
sex.means=cx.data.log %>%
group_by(SEX.mod) %>%
summarise(log.rCBR = mean(log.rCBR),
log.CX = mean(log.CX))
sex.means
cx.data.log$SEX.mod=factor(cx.data.log$SEX.mod,levels=
c("male","female"))
sex.CX.plot=ggplot(cx.data.log, aes(x=log.rCBR, y=log.CX, col=SEX.mod))+geom_point(alpha=0.3,size=4,stroke=NA)
sex.CX.plot
sex.CX.plot2=sex.CX.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
scale_color_manual(values = c("#59A3A2","#F0884D"))
sex.CX.plot2
sex.CX.plot3=sex.CX.plot2+geom_point(data=sex.means, size = 5)
sex.CX.plot3
# illustrates that there might be sex diffs, particulary when combined with phylogeny, but the biological meaningfulness is limited.
#see whether any substructures harbour sex effects more than others
#PB
sex.PB.plot=ggplot(cx.data, aes(x=log.rCBR, y=log.PB))+geom_point(aes(col=SEX.mod, size=3, alpha=0.5))
sex.PB.plot
sex.PB.plot2=sex.PB.plot+ theme_minimal()
sex.PB.plot2
#CBU
sex.CBU.plot=ggplot(cx.data, aes(x=log.rCBR, y=log.CBU))+geom_point(aes(col=SEX.mod, size=3, alpha=0.5))
sex.CBU.plot
sex.CBU.plot2=sex.CBU.plot+ theme_minimal()
sex.CBU.plot2
#CBL
sex.CBL.plot=ggplot(cx.data, aes(x=log.rCBR, y=log.CBL))+geom_point(aes(col=SEX.mod, size=3, alpha=0.5))
sex.CBL.plot
sex.CBL.plot2=sex.CBL.plot+ theme_minimal()
sex.CBL.plot2
#consistent sex effects here... can do posthoc tests but doesnt look too fascinating.
#we conclude that for the CX, we include sex as control factor but not examine biological context more closely.
#what about sex differences in terms of POTU and AOTU sizes?
#AOTU
#with new AOTU values
## define models
model.sex.aotu.A=log10(AOTU)~log10(rCBR)+SEX.mod
model.sex.aotu.B=log10(AOTU)~log10(rCBR)+SEX.mod+HelNonHel
model.sex.aotu.C=log10(AOTU)~log10(rCBR)+SEX.mod*HelNonHel
## run models
model.sex.aotu.A.1=MCMCglmm(model.sex.aotu.A, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.aotu.B.1=MCMCglmm(model.sex.aotu.B, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.aotu.C.1=MCMCglmm(model.sex.aotu.C, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the models twice
## run models
model.sex.aotu.A.2=MCMCglmm(model.sex.aotu.A, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.aotu.B.2=MCMCglmm(model.sex.aotu.B, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.aotu.C.2=MCMCglmm(model.sex.aotu.C, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
gelman.diag(mcmc.list(model.sex.aotu.A.1$Sol,model.sex.aotu.A.2$Sol))
gelman.diag(mcmc.list(model.sex.aotu.A.1$VCV,model.sex.aotu.A.2$VCV))
plot(mcmc.list(model.sex.aotu.A.1$VCV,model.sex.aotu.A.2$VCV))
autocorr(model.sex.aotu.A.1$Sol)
autocorr(model.sex.aotu.A.1$VCV)
gelman.diag(mcmc.list(model.sex.aotu.B.1$Sol,model.sex.aotu.B.2$Sol))
gelman.diag(mcmc.list(model.sex.aotu.B.1$VCV,model.sex.aotu.B.2$VCV))
plot(mcmc.list(model.sex.aotu.B.1$VCV,model.sex.aotu.B.2$VCV))
autocorr(model.sex.aotu.B.1$Sol)
autocorr(model.sex.aotu.B.1$VCV)
gelman.diag(mcmc.list(model.sex.aotu.C.1$Sol,model.sex.aotu.C.2$Sol))
gelman.diag(mcmc.list(model.sex.aotu.C.1$VCV,model.sex.aotu.C.2$VCV))
plot(mcmc.list(model.sex.aotu.C.1$VCV,model.sex.aotu.C.2$VCV))
autocorr(model.sex.aotu.C.1$Sol)
autocorr(model.sex.aotu.C.1$VCV)
DIC(model.sex.aotu.A.1)
DIC(model.sex.aotu.B.1)
DIC(model.sex.aotu.C.1)
#again simplest model.
summary(model.sex.aotu.A.1)
#highly signficiant.
sex.means.aotu=cx.data.log %>%
group_by(SEX.mod) %>%
summarise(log.rCBR = mean(log.rCBR),
log.AOTU = mean(log.AOTU))
sex.means.aotu
cx.data.log$SEX.mod=factor(cx.data.log$SEX.mod,levels=
c("male","female"))
sex.AOTU.plot=ggplot(cx.data.log, aes(x=log.rCBR, y=log.AOTU, col=SEX.mod))+geom_point(alpha=0.3,size=4,stroke=NA)
sex.AOTU.plot
sex.AOTU.plot2=sex.AOTU.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
scale_color_manual(values = c("#59A3A2","#F0884D"))
sex.AOTU.plot2
sex.AOTU.plot3=sex.AOTU.plot2+geom_point(data=sex.means.aotu, size =5)
sex.AOTU.plot3
#POTU
## define models
model.sex.potu.A=log10(POTU)~log10(rCBR)+SEX.mod
model.sex.potu.B=log10(POTU)~log10(rCBR)+SEX.mod+HelNonHel
model.sex.potu.C=log10(POTU)~log10(rCBR)+SEX.mod*HelNonHel
## run models
model.sex.potu.A.1=MCMCglmm(model.sex.potu.A, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.potu.B.1=MCMCglmm(model.sex.potu.B, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.potu.C.1=MCMCglmm(model.sex.potu.C, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the models twice
## run models
model.sex.potu.A.2=MCMCglmm(model.sex.potu.A, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.potu.B.2=MCMCglmm(model.sex.potu.B, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.sex.potu.C.2=MCMCglmm(model.sex.potu.C, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
gelman.diag(mcmc.list(model.sex.potu.A.1$Sol,model.sex.potu.A.2$Sol))
gelman.diag(mcmc.list(model.sex.potu.A.1$VCV,model.sex.potu.A.2$VCV))
plot(mcmc.list(model.sex.potu.A.1$VCV,model.sex.potu.A.2$VCV))
autocorr(model.sex.potu.A.1$Sol)
autocorr(model.sex.potu.A.1$VCV)
gelman.diag(mcmc.list(model.sex.potu.B.1$Sol,model.sex.potu.B.2$Sol))
gelman.diag(mcmc.list(model.sex.potu.B.1$VCV,model.sex.potu.B.2$VCV))
plot(mcmc.list(model.sex.potu.B.1$VCV,model.sex.potu.B.2$VCV))
autocorr(model.sex.potu.B.1$Sol)
autocorr(model.sex.potu.B.1$VCV)
gelman.diag(mcmc.list(model.sex.potu.C.1$Sol,model.sex.potu.C.2$Sol))
gelman.diag(mcmc.list(model.sex.potu.C.1$VCV,model.sex.potu.C.2$VCV))
plot(mcmc.list(model.sex.potu.C.1$VCV,model.sex.potu.C.2$VCV))
autocorr(model.sex.potu.C.1$Sol)
autocorr(model.sex.potu.C.1$VCV)
DIC(model.sex.potu.A.1)
DIC(model.sex.potu.B.1)
DIC(model.sex.potu.C.1)
#again simplest model.
summary(model.sex.potu.A.1)
#insignificant.
sex.means.POTU=cx.data.log %>%
group_by(SEX.mod) %>%
summarise(log.rCBR = mean(log.rCBR),
log.POTU = mean(log.POTU))
sex.means.POTU
cx.data.log$SEX.mod=factor(cx.data.log$SEX.mod,levels=
c("male","female"))
sex.POTU.plot=ggplot(cx.data.log, aes(x=log.rCBR, y=log.POTU, col=SEX.mod))+geom_point(alpha=0.3,size=4,stroke=NA)
sex.POTU.plot
sex.POTU.plot2=sex.POTU.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),panel.background = element_blank(), axis.line = element_line(colour = "black"))+
scale_color_manual(values = c("#59A3A2","#F0884D"))
sex.POTU.plot2
sex.POTU.plot3=sex.POTU.plot2+geom_point(data=sex.means.POTU, size = 5)
sex.POTU.plot3
sex.plot=ggarrange(sex.CX.plot3,sex.AOTU.plot3,sex.POTU.plot3,
ncol = 3, nrow = 1,common.legend = TRUE, legend="bottom")
sex.plot
##### CLADE DIFFERENCES (related to Figure 2 and S2) #####
#number of individuals per species.
count.data=cx.data %>% count(phylo)
#TREE
tree = read.nexus("Heliconiini.trees")
tree=force.ultrametric(tree,method="nnls")
plot(tree)
inv.phylo=inverseA(tree,nodes="TIPS",scale=TRUE)
##MCMCglmm
# set priors
prior=list(G=list(G1=list(V=1,nu=0.02)),R=list(V=1,nu=0.02)) #standard priors as explored by SHM
## define model
model.CX.total=log10(Total.CX.woNO)~log10(rCBR)+PollenF+SEX.mod
## run model
model.CX.1=MCMCglmm(model.CX.total, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CX.2=MCMCglmm(model.CX.total, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CX.1$Sol,model.CX.2$Sol))
gelman.diag(mcmc.list(model.CX.1$VCV,model.CX.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CX.1$VCV,model.CX.2$VCV))
plot(mcmc.list(model.CX.1$Sol,model.CX.2$Sol))
#all good.
summary(model.CX.1)
#non significant.
model.AOTU=log10(AOTU)~log10(rCBR)+PollenF+SEX.mod
model.AOTU.1=MCMCglmm(model.AOTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.AOTU.2=MCMCglmm(model.AOTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.AOTU.1$Sol,model.AOTU.2$Sol))
gelman.diag(mcmc.list(model.AOTU.1$VCV,model.AOTU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.AOTU.1$VCV,model.AOTU.2$VCV))
plot(mcmc.list(model.AOTU.1$Sol,model.AOTU.2$Sol))
#all good.
summary(model.AOTU.1)
#insignificant.
model.POTU=log10(POTU)~log10(rCBR)+PollenF+SEX.mod
model.POTU.1=MCMCglmm(model.POTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.POTU.2=MCMCglmm(model.POTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.POTU.1$Sol,model.POTU.2$Sol))
gelman.diag(mcmc.list(model.POTU.1$VCV,model.POTU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.POTU.1$VCV,model.POTU.2$VCV))
plot(mcmc.list(model.POTU.1$Sol,model.POTU.2$Sol))
#all good.
summary(model.POTU.1)
#insignificant.
model.PB=log10(PB)~log10(rCBR)+PollenF+SEX.mod
model.PB.1=MCMCglmm(model.PB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.PB.2=MCMCglmm(model.PB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.PB.1$Sol,model.PB.2$Sol))
gelman.diag(mcmc.list(model.PB.1$VCV,model.PB.2$VCV))
#random and fixed okay
plot(mcmc.list(model.PB.1$VCV,model.PB.2$VCV))
plot(mcmc.list(model.PB.1$Sol,model.PB.2$Sol))
#all good.
summary(model.PB.1)
#insignificant.
model.CBL=log10(CBL)~log10(rCBR)+PollenF+SEX.mod
model.CBL.1=MCMCglmm(model.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBL.2=MCMCglmm(model.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBL.1$Sol,model.CBL.2$Sol))
gelman.diag(mcmc.list(model.CBL.1$VCV,model.CBL.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.1$VCV,model.CBL.2$VCV))
plot(mcmc.list(model.CBL.1$Sol,model.CBL.2$Sol))
#all good.
summary(model.CBL.1)
#insignificant.
model.CBU=log10(CBU)~log10(rCBR)+PollenF+SEX.mod
model.CBU.1=MCMCglmm(model.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBU.2=MCMCglmm(model.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBU.1$Sol,model.CBU.2$Sol))
gelman.diag(mcmc.list(model.CBU.1$VCV,model.CBU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBU.1$VCV,model.CBU.2$VCV))
plot(mcmc.list(model.CBU.1$Sol,model.CBU.2$Sol))
#all good.
summary(model.CBU.1)
#insignificant.
model.NO=log10(NO)~log10(rCBR)+PollenF+SEX.mod
model.NO.1=MCMCglmm(model.NO, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.2=MCMCglmm(model.NO, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.1$Sol,model.NO.2$Sol))
gelman.diag(mcmc.list(model.NO.1$VCV,model.NO.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.1$VCV,model.NO.2$VCV))
plot(mcmc.list(model.NO.1$Sol,model.NO.2$Sol))
#all good.
summary(model.NO.1)
#insignificant.
CX.plot=ggplot(cx.data, aes(x=log10(rCBR), y=log10(Total.CX.woNO),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
CX.plot
CX.plot2=CX.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
CX.plot2
#ggplotly(CX.plot2)
#speciesaverages.
spec.avg=read.table(file="SpecAvg-plot_vs2.txt", sep="\t", header=TRUE)
str(spec.avg)
CX.plot3=CX.plot2+geom_point(data=spec.avg,size=5)
CX.plot3
AOTU.plot=ggplot(cx.data, aes(x=log10(rCBR), y=log10(AOTU),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
AOTU.plot
AOTU.plot2=AOTU.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
AOTU.plot2
AOTU.plot3=AOTU.plot2+geom_point(data=spec.avg,size=5)
AOTU.plot3
POTU.plot=ggplot(cx.data, aes(x=log10(rCBR), y=log10(POTU),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
POTU.plot
POTU.plot2=POTU.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
POTU.plot2
POTU.plot3=POTU.plot2+geom_point(data=spec.avg,size=5)
POTU.plot3
CBU.plot=ggplot(cx.data, aes(x=log10(rCBR), y=log10(CBU),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
CBU.plot
CBU.plot2=CBU.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
CBU.plot2
CBU.plot3=CBU.plot2+geom_point(data=spec.avg,size=5)
CBU.plot3
CBL.plot=ggplot(cx.data, aes(x=log10(rCBR), y=log10(CBL),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
CBL.plot
CBL.plot2=CBL.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
CBL.plot2
CBL.plot3=CBL.plot2+geom_point(data=spec.avg,size=5)
CBL.plot3
PB.plot=ggplot(cx.data, aes(x=log10(rCBR), y=log10(PB),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
PB.plot
PB.plot2=PB.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
PB.plot2
PB.plot3=PB.plot2+geom_point(data=spec.avg,size=5)
PB.plot3
NO.data=subset(cx.data, NO >= 20)
NO.data.avg=subset(spec.avg,NO >= 20)
NO.plot=ggplot(NO.data, aes(x=log10(rCBR), y=log10(NO),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
NO.plot
NO.plot2=NO.plot+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
NO.plot2
NO.plot3=NO.plot2+geom_point(data=NO.data.avg,size=5)
NO.plot3
pollen.p1= ggarrange(PB.plot3,CBU.plot3, CBL.plot3 ,NO.plot3,
ncol = 4, nrow = 1,common.legend = TRUE, legend="bottom" )
pollen.p1
pollen.p2=ggarrange(CX.plot3, AOTU.plot3, POTU.plot3,
ncol = 3, nrow = 1,common.legend = TRUE, legend="bottom" )
pollen.p2
xplot=ggdensity(cx.data, "log10(rCBR)", fill = "Grp.color")+scale_fill_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))+clean_theme()
xplot
CX.yplot=ggdensity(cx.data, "log10(Total.CX.woNO)", fill = "Grp.color")+scale_fill_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))+coord_flip()+clean_theme()
CX.yplot
CX.plot.full=ggarrange(xplot, NULL, CX.plot3, CX.yplot,
ncol = 2, nrow = 2, align = "hv",
widths = c(2, 1), heights = c(1, 2),
common.legend = TRUE)
CX.plot.full
##### MB-CX co-evolution DIFFERENCES (related to Figure 2 and S2) #####
#TREE
tree = read.nexus("Heliconiini.trees")
tree=force.ultrametric(tree,method="nnls")
plot(tree)
inv.phylo=inverseA(tree,nodes="TIPS",scale=TRUE)
##MCMCglmm
# set priors
prior=list(G=list(G1=list(V=1,nu=0.02)),R=list(V=1,nu=0.02)) #standard priors as explored by SHM
#total CX
## define model
model.CX.MB.PF=log10(Total.CX.woNO)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.CX.MB.PF.1=MCMCglmm(model.CX.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CX.MB.PF.2=MCMCglmm(model.CX.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CX.MB.PF.1$Sol,model.CX.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.CX.MB.PF.1$VCV,model.CX.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CX.MB.PF.1$VCV,model.CX.MB.PF.2$VCV))
plot(mcmc.list(model.CX.MB.PF.1$Sol,model.CX.MB.PF.2$Sol))
#all good.
summary(model.CX.MB.PF.1)
#MB and rCBR significant.
model.CX.MB.PFn=log10(Total.CX.woNO)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## run model
model.CX.MB.PFn.1=MCMCglmm(model.CX.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CX.MB.PFn.2=MCMCglmm(model.CX.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CX.MB.PFn.1$Sol,model.CX.MB.PFn.2$Sol))
gelman.diag(mcmc.list(model.CX.MB.PFn.1$VCV,model.CX.MB.PFn.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CX.MB.PFn.1$VCV,model.CX.MB.PFn.2$VCV))
plot(mcmc.list(model.CX.MB.PFn.1$Sol,model.CX.MB.PFn.2$Sol))
#all good.
DIC(model.CX.MB.PFn.1)
DIC(model.CX.MB.PF.1)
#very similar.
sum.CX=summary(model.CX.MB.PF.1)
sum.CX.output=sum.CX$solutions
#total CX. clade effects.
## define model
model.CX.MB.cl=log10(Total.CX.woNO)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## run model
model.CX.MB.cl.1=MCMCglmm(model.CX.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CX.MB.cl.2=MCMCglmm(model.CX.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CX.MB.cl.1$Sol,model.CX.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.CX.MB.cl.1$VCV,model.CX.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CX.MB.cl.1$VCV,model.CX.MB.cl.2$VCV))
plot(mcmc.list(model.CX.MB.cl.1$Sol,model.CX.MB.cl.2$Sol))
#all good.
summary(model.CX.MB.cl.1)
#MB and rCBR significant.
#comparison to null model.
model.CX.MB.cln=log10(Total.CX.woNO)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## run model
model.CX.MB.cln.1=MCMCglmm(model.CX.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CX.MB.cln.2=MCMCglmm(model.CX.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CX.MB.cln.1$Sol,model.CX.MB.cln.2$Sol))
gelman.diag(mcmc.list(model.CX.MB.cln.1$VCV,model.CX.MB.cln.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CX.MB.cln.1$VCV,model.CX.MB.cln.2$VCV))
plot(mcmc.list(model.CX.MB.cln.1$Sol,model.CX.MB.cln.2$Sol))
DIC(model.CX.MB.cln.1)
DIC(model.CX.MB.cl.1)
#AOTU
## define model
model.AOTU.MB.PF=log10(AOTU)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.AOTU.MB.PF.1=MCMCglmm(model.AOTU.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.AOTU.MB.PF.2=MCMCglmm(model.AOTU.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.AOTU.MB.PF.1$Sol,model.AOTU.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.AOTU.MB.PF.1$VCV,model.AOTU.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.AOTU.MB.PF.1$VCV,model.AOTU.MB.PF.2$VCV))
plot(mcmc.list(model.AOTU.MB.PF.1$Sol,model.AOTU.MB.PF.2$Sol))
#all good.
## null model
model.AOTU.MB.PFn=log10(AOTU)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## model
model.AOTU.MB.PFn.1=MCMCglmm(model.AOTU.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.AOTU.MB.PFn.2=MCMCglmm(model.AOTU.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.AOTU.MB.PF.1)
DIC(model.AOTU.MB.PFn.1)
summary(model.AOTU.MB.PF.1)
#MB and rCBR significant.
sum.aotu=summary(model.AOTU.MB.PF.1)
sum.aotu.output=sum.aotu$solutions
#clade diffs
model.AOTU.MB.cl=log10(AOTU)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## model
model.AOTU.MB.cl.1=MCMCglmm(model.AOTU.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.AOTU.MB.cl.2=MCMCglmm(model.AOTU.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.AOTU.MB.cl.1$Sol,model.AOTU.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.AOTU.MB.cl.1$VCV,model.AOTU.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.AOTU.MB.cl.1$VCV,model.AOTU.MB.cl.2$VCV))
plot(mcmc.list(model.AOTU.MB.cl.1$Sol,model.AOTU.MB.cl.2$Sol))
#all good.
#null model
model.AOTU.MB.cln=log10(AOTU)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## model
model.AOTU.MB.cln.1=MCMCglmm(model.AOTU.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.AOTU.MB.cln.2=MCMCglmm(model.AOTU.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.AOTU.MB.cl.1)
DIC(model.AOTU.MB.cln.1)
summary(model.AOTU.MB.cl.1)
#MB and rCBR significant.
sum.aotu.cl=summary(model.AOTU.MB.cl.1)
sum.aotu.cl.output=sum.aotu.cl$solutions
#POTU
## define model
model.POTU.MB.PF=log10(POTU)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.POTU.MB.PF.1=MCMCglmm(model.POTU.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.POTU.MB.PF.2=MCMCglmm(model.POTU.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.POTU.MB.PF.1$Sol,model.POTU.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.POTU.MB.PF.1$VCV,model.POTU.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.POTU.MB.PF.1$VCV,model.POTU.MB.PF.2$VCV))
plot(mcmc.list(model.POTU.MB.PF.1$Sol,model.POTU.MB.PF.2$Sol))
#all good.
## null model
model.POTU.MB.PFn=log10(POTU)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## model
model.POTU.MB.PFn.1=MCMCglmm(model.POTU.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.POTU.MB.PFn.2=MCMCglmm(model.POTU.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.POTU.MB.PF.1)
DIC(model.POTU.MB.PFn.1)
summary(model.POTU.MB.PF.1)
#MB and rCBR significant.
sum.POTU=summary(model.POTU.MB.PF.1)
sum.POTU.output=sum.POTU$solutions
#clade diffs
model.POTU.MB.cl=log10(POTU)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## model
model.POTU.MB.cl.1=MCMCglmm(model.POTU.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.POTU.MB.cl.2=MCMCglmm(model.POTU.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.POTU.MB.cl.1$Sol,model.POTU.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.POTU.MB.cl.1$VCV,model.POTU.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.POTU.MB.cl.1$VCV,model.POTU.MB.cl.2$VCV))
plot(mcmc.list(model.POTU.MB.cl.1$Sol,model.POTU.MB.cl.2$Sol))
#all good.
#null model
model.POTU.MB.cln=log10(POTU)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## model
model.POTU.MB.cln.1=MCMCglmm(model.POTU.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.POTU.MB.cln.2=MCMCglmm(model.POTU.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.POTU.MB.cl.1)
DIC(model.POTU.MB.cln.1)
summary(model.POTU.MB.cl.1)
#MB and rCBR significant.
sum.POTU.cl=summary(model.POTU.MB.cl.1)
sum.POTU.cl.output=sum.POTU.cl$solutions
#CBU
## define model
model.CBU.MB.PF=log10(CBU)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.CBU.MB.PF.1=MCMCglmm(model.CBU.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBU.MB.PF.2=MCMCglmm(model.CBU.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBU.MB.PF.1$Sol,model.CBU.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.CBU.MB.PF.1$VCV,model.CBU.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBU.MB.PF.1$VCV,model.CBU.MB.PF.2$VCV))
plot(mcmc.list(model.CBU.MB.PF.1$Sol,model.CBU.MB.PF.2$Sol))
#all good.
## null model
model.CBU.MB.PFn=log10(CBU)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## model
model.CBU.MB.PFn.1=MCMCglmm(model.CBU.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBU.MB.PFn.2=MCMCglmm(model.CBU.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.CBU.MB.PF.1)
DIC(model.CBU.MB.PFn.1)
summary(model.CBU.MB.PF.1)
#MB and rCBR significant.
sum.CBU=summary(model.CBU.MB.PF.1)
sum.CBU.output=sum.CBU$solutions
#clade diffs
model.CBU.MB.cl=log10(CBU)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## model
model.CBU.MB.cl.1=MCMCglmm(model.CBU.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBU.MB.cl.2=MCMCglmm(model.CBU.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBU.MB.cl.1$Sol,model.CBU.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.CBU.MB.cl.1$VCV,model.CBU.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBU.MB.cl.1$VCV,model.CBU.MB.cl.2$VCV))
plot(mcmc.list(model.CBU.MB.cl.1$Sol,model.CBU.MB.cl.2$Sol))
#all good.
#null model
model.CBU.MB.cln=log10(CBU)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## model
model.CBU.MB.cln.1=MCMCglmm(model.CBU.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBU.MB.cln.2=MCMCglmm(model.CBU.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.CBU.MB.cl.1)
DIC(model.CBU.MB.cln.1)
summary(model.CBU.MB.cl.1)
#MB and rCBR significant.
sum.CBU.cl=summary(model.CBU.MB.cl.1)
sum.CBU.cl.output=sum.CBU.cl$solutions
#CBL
## define model
model.CBL.MB.PF=log10(CBL)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.CBL.MB.PF.1=MCMCglmm(model.CBL.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBL.MB.PF.2=MCMCglmm(model.CBL.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBL.MB.PF.1$Sol,model.CBL.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.CBL.MB.PF.1$VCV,model.CBL.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.MB.PF.1$VCV,model.CBL.MB.PF.2$VCV))
plot(mcmc.list(model.CBL.MB.PF.1$Sol,model.CBL.MB.PF.2$Sol))
#all good.
## null model
model.CBL.MB.PFn=log10(CBL)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## model
model.CBL.MB.PFn.1=MCMCglmm(model.CBL.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBL.MB.PFn.2=MCMCglmm(model.CBL.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.CBL.MB.PF.1)
DIC(model.CBL.MB.PFn.1)
summary(model.CBL.MB.PF.1)
#MB and rCBR significant.
sum.CBL=summary(model.CBL.MB.PF.1)
sum.CBL.output=sum.CBL$solutions
#clade diffs
model.CBL.MB.cl=log10(CBL)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## model
model.CBL.MB.cl.1=MCMCglmm(model.CBL.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBL.MB.cl.2=MCMCglmm(model.CBL.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBL.MB.cl.1$Sol,model.CBL.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.CBL.MB.cl.1$VCV,model.CBL.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.MB.cl.1$VCV,model.CBL.MB.cl.2$VCV))
plot(mcmc.list(model.CBL.MB.cl.1$Sol,model.CBL.MB.cl.2$Sol))
#all good.
#null model
model.CBL.MB.cln=log10(CBL)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## model
model.CBL.MB.cln.1=MCMCglmm(model.CBL.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBL.MB.cln.2=MCMCglmm(model.CBL.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.CBL.MB.cl.1)
DIC(model.CBL.MB.cln.1)
summary(model.CBL.MB.cl.1)
#MB and rCBR significant.
sum.CBL.cl=summary(model.CBL.MB.cl.1)
sum.CBL.cl.output=sum.CBL.cl$solutions
#PB
## define model
model.PB.MB.PF=log10(PB)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.PB.MB.PF.1=MCMCglmm(model.PB.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.PB.MB.PF.2=MCMCglmm(model.PB.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.PB.MB.PF.1$Sol,model.PB.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.PB.MB.PF.1$VCV,model.PB.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.PB.MB.PF.1$VCV,model.PB.MB.PF.2$VCV))
plot(mcmc.list(model.PB.MB.PF.1$Sol,model.PB.MB.PF.2$Sol))
#all good.
## null model
model.PB.MB.PFn=log10(PB)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## model
model.PB.MB.PFn.1=MCMCglmm(model.PB.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.PB.MB.PFn.2=MCMCglmm(model.PB.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.PB.MB.PF.1)
DIC(model.PB.MB.PFn.1)
summary(model.PB.MB.PF.1)
#MB and rCBR significant.
sum.PB=summary(model.PB.MB.PF.1)
sum.PB.output=sum.PB$solutions
#clade diffs
model.PB.MB.cl=log10(PB)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## model
model.PB.MB.cl.1=MCMCglmm(model.PB.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.PB.MB.cl.2=MCMCglmm(model.PB.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.PB.MB.cl.1$Sol,model.PB.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.PB.MB.cl.1$VCV,model.PB.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.PB.MB.cl.1$VCV,model.PB.MB.cl.2$VCV))
plot(mcmc.list(model.PB.MB.cl.1$Sol,model.PB.MB.cl.2$Sol))
#all good.
#null model
model.PB.MB.cln=log10(PB)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## model
model.PB.MB.cln.1=MCMCglmm(model.PB.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.PB.MB.cln.2=MCMCglmm(model.PB.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.PB.MB.cl.1)
DIC(model.PB.MB.cln.1)
summary(model.PB.MB.cl.1)
#MB and rCBR significant.
sum.PB.cl=summary(model.PB.MB.cl.1)
sum.PB.cl.output=sum.PB.cl$solutions
#NO
## define model
model.NO.MB.PF=log10(NO)~log10(rCBR)+log10(total.MB)*PollenF+SEX.mod
## run model
model.NO.MB.PF.1=MCMCglmm(model.NO.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.MB.PF.2=MCMCglmm(model.NO.MB.PF, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.MB.PF.1$Sol,model.NO.MB.PF.2$Sol))
gelman.diag(mcmc.list(model.NO.MB.PF.1$VCV,model.NO.MB.PF.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.MB.PF.1$VCV,model.NO.MB.PF.2$VCV))
plot(mcmc.list(model.NO.MB.PF.1$Sol,model.NO.MB.PF.2$Sol))
#all good.
## null model
model.NO.MB.PFn=log10(NO)~log10(rCBR)+log10(total.MB)+PollenF+SEX.mod
## model
model.NO.MB.PFn.1=MCMCglmm(model.NO.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.MB.PFn.2=MCMCglmm(model.NO.MB.PFn, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.NO.MB.PF.1)
DIC(model.NO.MB.PFn.1)
summary(model.NO.MB.PF.1)
#MB and rCBR significant.
sum.NO=summary(model.NO.MB.PF.1)
sum.NO.output=sum.NO$solutions
#clade diffs
model.NO.MB.cl=log10(NO)~log10(rCBR)+log10(total.MB)*Grp.color+SEX.mod
## model
model.NO.MB.cl.1=MCMCglmm(model.NO.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.MB.cl.2=MCMCglmm(model.NO.MB.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.MB.cl.1$Sol,model.NO.MB.cl.2$Sol))
gelman.diag(mcmc.list(model.NO.MB.cl.1$VCV,model.NO.MB.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.MB.cl.1$VCV,model.NO.MB.cl.2$VCV))
plot(mcmc.list(model.NO.MB.cl.1$Sol,model.NO.MB.cl.2$Sol))
#all good.
#null model
model.NO.MB.cln=log10(NO)~log10(rCBR)+log10(total.MB)+Grp.color+SEX.mod
## model
model.NO.MB.cln.1=MCMCglmm(model.NO.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.MB.cln.2=MCMCglmm(model.NO.MB.cln, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
DIC(model.NO.MB.cl.1)
DIC(model.NO.MB.cln.1)
summary(model.NO.MB.cl.1)
#MB and rCBR significant.
sum.NO.cl=summary(model.NO.MB.cl.1)
sum.NO.cl.output=sum.NO.cl$solutions
#plotting
CX.plot.MB=ggplot(cx.data, aes(x=log10(total.MB), y=log10(Total.CX.woNO),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
CX.plot.MB
CX.plot.MB.2=CX.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
CX.plot.MB.2
#ggplotly(CX.plot2)
#speciesaverages.
spec.avg=read.table(file="SpecAvg-plot.txt", sep="\t", header=TRUE)
str(spec.avg)
CX.plot.MB.3=CX.plot.MB.2+geom_point(data=spec.avg,size=5)
CX.plot.MB.3
AOTU.plot.MB=ggplot(cx.data, aes(x=log10(total.MB), y=log10(AOTU),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
AOTU.plot.MB
AOTU.plot.MB.2=AOTU.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
AOTU.plot.MB.2
AOTU.plot.MB.3=AOTU.plot.MB.2+geom_point(data=spec.avg,size=5)
AOTU.plot.MB.3
POTU.plot.MB=ggplot(cx.data, aes(x=log10(total.MB), y=log10(POTU),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
POTU.plot.MB
POTU.plot.MB.2=POTU.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
POTU.plot.MB.2
POTU.plot.MB.3=POTU.plot.MB.2+geom_point(data=spec.avg,size=5)
POTU.plot.MB.3
CBU.plot.MB=ggplot(cx.data, aes(x=log10(total.MB), y=log10(CBU),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
CBU.plot.MB
CBU.plot.MB.2=CBU.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
CBU.plot.MB.2
CBU.plot.MB.3=CBU.plot.MB.2+geom_point(data=spec.avg,size=5)
CBU.plot.MB.3
CBL.plot.MB=ggplot(cx.data, aes(x=log10(total.MB), y=log10(CBL),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
CBL.plot.MB
CBL.plot.MB.2=CBL.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
CBL.plot.MB.2
CBL.plot.MB.3=CBL.plot.MB.2+geom_point(data=spec.avg,size=5)
CBL.plot.MB.3
PB.plot.MB=ggplot(cx.data, aes(x=log10(total.MB), y=log10(PB),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
PB.plot.MB
PB.plot.MB.2=PB.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
PB.plot.MB.2
PB.plot.MB.3=PB.plot.MB.2+geom_point(data=spec.avg,size=5)
PB.plot.MB.3
NO.data=subset(cx.data, NO >= 20)
NO.data.avg=subset(spec.avg,NO >= 20)
NO.plot.MB=ggplot(NO.data, aes(x=log10(total.MB), y=log10(NO),col=Grp.color,size=2.5))+geom_point(alpha=0.3,stroke=NA)
NO.plot.MB
NO.plot.MB.2=NO.plot.MB+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))+scale_colour_manual(values = c("#7d9c2c","#84478e", "#c64f37","#ed6d26","#1a90b3"))
NO.plot.MB.2
NO.plot.MB.3=NO.plot.MB.2+geom_point(data=NO.data.avg,size=5)
NO.plot.MB.3
pollen.MB.p1= ggarrange(PB.plot.MB.3,CBU.plot.MB.3, CBL.plot.MB.3 ,NO.plot.MB.3,
ncol = 4, nrow = 1,common.legend = TRUE, legend="bottom" )
pollen.MB.p1
pollen.MB.p2=ggarrange(CX.plot.MB.3, AOTU.plot.MB.3, POTU.plot.MB.3,
ncol = 3, nrow = 1,common.legend = TRUE, legend="bottom" )
pollen.MB.p2
#single pair models.
#TREE
tree = read.nexus("Heliconiini.trees")
tree=force.ultrametric(tree,method="nnls")
plot(tree)
inv.phylo=inverseA(tree,nodes="TIPS",scale=TRUE)
##MCMCglmm
# set priors
prior=list(G=list(G1=list(V=1,nu=0.02)),R=list(V=1,nu=0.02)) #standard priors as explored by SHM
#CBU tNO
## define model
model.NO.CBU=log10(NO)~log10(rCBR)+log10(CBU)*PollenF+SEX.mod
## run model
model.NO.CBU.1=MCMCglmm(model.NO.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.CBU.2=MCMCglmm(model.NO.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.CBU.1$Sol,model.NO.CBU.2$Sol))
gelman.diag(mcmc.list(model.NO.CBU.1$VCV,model.NO.CBU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.CBU.1$VCV,model.NO.CBU.2$VCV))
plot(mcmc.list(model.NO.CBU.1$Sol,model.NO.CBU.2$Sol))
#all good.
summary(model.NO.CBU.1)
#hm no diffs according to PF.
#CBU tNO
## define model
model.NO.CBU.cl=log10(NO)~log10(rCBR)+log10(CBU)*Grp.color+SEX.mod
## run model
model.NO.CBU.cl.1=MCMCglmm(model.NO.CBU.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.CBU.cl.2=MCMCglmm(model.NO.CBU.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.CBU.cl.1$Sol,model.NO.CBU.cl.2$Sol))
gelman.diag(mcmc.list(model.NO.CBU.cl.1$VCV,model.NO.CBU.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.CBU.cl.1$VCV,model.NO.CBU.cl.2$VCV))
plot(mcmc.list(model.NO.CBU.cl.1$Sol,model.NO.CBU.cl.2$Sol))
#all good.
summary(model.NO.CBU.cl.1)
#CBL tNO
## define model
model.NO.CBL=log10(NO)~log10(rCBR)+log10(CBL)*PollenF+SEX.mod
## run model
model.NO.CBL.1=MCMCglmm(model.NO.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.CBL.2=MCMCglmm(model.NO.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.CBL.1$Sol,model.NO.CBL.2$Sol))
gelman.diag(mcmc.list(model.NO.CBL.1$VCV,model.NO.CBL.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.CBL.1$VCV,model.NO.CBL.2$VCV))
plot(mcmc.list(model.NO.CBL.1$Sol,model.NO.CBL.2$Sol))
#all good.
summary(model.NO.CBL.1)
#hm no diffs according to PF.
#CBL tNO
## define model
model.NO.CBL.cl=log10(NO)~log10(rCBR)+log10(CBL)*Grp.color+SEX.mod
## run model
model.NO.CBL.cl.1=MCMCglmm(model.NO.CBL.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.CBL.cl.2=MCMCglmm(model.NO.CBL.cl, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.CBL.cl.1$Sol,model.NO.CBL.cl.2$Sol))
gelman.diag(mcmc.list(model.NO.CBL.cl.1$VCV,model.NO.CBL.cl.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.CBL.cl.1$VCV,model.NO.CBL.cl.2$VCV))
plot(mcmc.list(model.NO.CBL.cl.1$Sol,model.NO.CBL.cl.2$Sol))
#all good.
summary(model.NO.CBL.cl.1)
##### CX, AOTU, POTU within DIFFERENCES (Figure 3) #####
setwd("C:/R-Analysis")
cx.data=read.table(file="Hel-CX-MB_vs2.txt", sep="\t", header=TRUE)
str(cx.data)
#TREE
tree = read.nexus("Heliconiini.trees")
tree=force.ultrametric(tree,method="nnls")
plot(tree)
inv.phylo=inverseA(tree,nodes="TIPS",scale=TRUE)
##MCMCglmm
# set priors
prior=list(G=list(G1=list(V=1,nu=0.02)),R=list(V=1,nu=0.02)) #standard priors as explored by SHM
#AOTU to rest
## define model
model.AOTU.win=log10(AOTU)~log10(rCBR)+log10(POTU)+log10(CBU)+log10(CBL)+log10(PB)+SEX.mod
## run model
model.AOTU.win.1=MCMCglmm(model.AOTU.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.AOTU.win.2=MCMCglmm(model.AOTU.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.AOTU.win.1$Sol,model.AOTU.win.2$Sol))
gelman.diag(mcmc.list(model.AOTU.win.1$VCV,model.AOTU.win.2$VCV))
#random and fixed okay
plot(mcmc.list(model.AOTU.win.1$VCV,model.AOTU.win.2$VCV))
plot(mcmc.list(model.AOTU.win.1$Sol,model.AOTU.win.2$Sol))
#all good.
sum.AOTU=summary(model.AOTU.win.1)
sum.AOTU.output=sum.AOTU$solutions
#POTU to rest
## define model
model.POTU.win=log10(POTU)~log10(rCBR)+log10(AOTU)+log10(CBU)+log10(CBL)+log10(PB)+SEX.mod
## run model
model.POTU.win.1=MCMCglmm(model.POTU.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.POTU.win.2=MCMCglmm(model.POTU.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.POTU.win.1$Sol,model.POTU.win.2$Sol))
gelman.diag(mcmc.list(model.POTU.win.1$VCV,model.POTU.win.2$VCV))
#random and fixed okay
plot(mcmc.list(model.POTU.win.1$VCV,model.POTU.win.2$VCV))
plot(mcmc.list(model.POTU.win.1$Sol,model.POTU.win.2$Sol))
#all good.
sum.POTU=summary(model.POTU.win.1)
sum.POTU.output=sum.POTU$solutions
#CBU to rest
## define model
model.CBU.win=log10(CBU)~log10(rCBR)+log10(AOTU)+log10(POTU)+log10(CBL)+log10(PB)+SEX.mod
## run model
model.CBU.win.1=MCMCglmm(model.CBU.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBU.win.2=MCMCglmm(model.CBU.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBU.win.1$Sol,model.CBU.win.2$Sol))
gelman.diag(mcmc.list(model.CBU.win.1$VCV,model.CBU.win.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBU.win.1$VCV,model.CBU.win.2$VCV))
plot(mcmc.list(model.CBU.win.1$Sol,model.CBU.win.2$Sol))
#all good.
sum.CBU=summary(model.CBU.win.1)
sum.CBU.output=sum.CBU$solutions
#CBL to rest
## define model
model.CBL.win=log10(CBL)~log10(rCBR)+log10(AOTU)+log10(POTU)+log10(CBU)+log10(PB)+SEX.mod
## run model
model.CBL.win.1=MCMCglmm(model.CBL.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.CBL.win.2=MCMCglmm(model.CBL.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.CBL.win.1$Sol,model.CBL.win.2$Sol))
gelman.diag(mcmc.list(model.CBL.win.1$VCV,model.CBL.win.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.win.1$VCV,model.CBL.win.2$VCV))
plot(mcmc.list(model.CBL.win.1$Sol,model.CBL.win.2$Sol))
#all good.
sum.CBL=summary(model.CBL.win.1)
sum.CBL.output=sum.CBL$solutions
#PB to rest
## define model
model.PB.win=log10(PB)~log10(rCBR)+log10(AOTU)+log10(POTU)+log10(CBU)+log10(CBL)+SEX.mod
## run model
model.PB.win.1=MCMCglmm(model.PB.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.PB.win.2=MCMCglmm(model.PB.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.PB.win.1$Sol,model.PB.win.2$Sol))
gelman.diag(mcmc.list(model.PB.win.1$VCV,model.PB.win.2$VCV))
#random and fixed okay
plot(mcmc.list(model.PB.win.1$VCV,model.PB.win.2$VCV))
plot(mcmc.list(model.PB.win.1$Sol,model.PB.win.2$Sol))
#all good.
sum.PB=summary(model.PB.win.1)
sum.PB.output=sum.PB$solutions
#NO to rest
NO.data=subset(cx.data, NO >= 20)
## define model
model.NO.win=log10(NO)~log10(rCBR)+log10(AOTU)+log10(POTU)+log10(CBU)+log10(CBL)+log10(PB)+SEX.mod
## run model
model.NO.win.1=MCMCglmm(model.NO.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
## run the model twice
model.NO.win.2=MCMCglmm(model.NO.win, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=NO.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.NO.win.1$Sol,model.NO.win.2$Sol))
gelman.diag(mcmc.list(model.NO.win.1$VCV,model.NO.win.2$VCV))
#random and fixed okay
plot(mcmc.list(model.NO.win.1$VCV,model.NO.win.2$VCV))
plot(mcmc.list(model.NO.win.1$Sol,model.NO.win.2$Sol))
#all good.
sum.NO=summary(model.NO.win.1)
sum.NO.output=sum.NO$solutions
#rerun significant relationships to test effects of pollen feeding
#model
model.AOTU.CBU=log10(AOTU)~log10(rCBR)+log10(CBU)*PollenF+SEX.mod
model.AOTU.CBL=log10(AOTU)~log10(rCBR)+log10(CBL)*PollenF+SEX.mod
model.POTU.PB=log10(POTU)~log10(rCBR)+log10(PB)*PollenF+SEX.mod
model.CBU.AOTU=log10(CBU)~log10(rCBR)+log10(AOTU)*PollenF+SEX.mod
model.CBU.CBL=log10(CBU)~log10(rCBR)+log10(CBL)*PollenF+SEX.mod
model.CBL.AOTU=log10(CBL)~log10(rCBR)+log10(AOTU)*PollenF+SEX.mod
model.CBL.CBU=log10(CBL)~log10(rCBR)+log10(CBU)*PollenF+SEX.mod
model.CBL.PB=log10(CBL)~log10(rCBR)+log10(PB)*PollenF+SEX.mod
model.PB.POTU=log10(PB)~log10(rCBR)+log10(POTU)*PollenF+SEX.mod
model.PB.CBL=log10(PB)~log10(rCBR)+log10(CBL)*PollenF+SEX.mod
## run model
model.AOTU.CBU.1=MCMCglmm(model.AOTU.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.AOTU.CBL.1=MCMCglmm(model.AOTU.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.POTU.PB.1=MCMCglmm(model.POTU.PB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBU.AOTU.1=MCMCglmm(model.CBU.AOTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBU.CBL.1=MCMCglmm(model.CBU.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBL.AOTU.1=MCMCglmm(model.CBL.AOTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBL.CBU.1=MCMCglmm(model.CBL.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBL.PB.1=MCMCglmm(model.CBL.PB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.PB.POTU.1=MCMCglmm(model.PB.POTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.PB.CBL.1=MCMCglmm(model.PB.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
## run model again
model.AOTU.CBU.2=MCMCglmm(model.AOTU.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.AOTU.CBL.2=MCMCglmm(model.AOTU.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.POTU.PB.2=MCMCglmm(model.POTU.PB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBU.AOTU.2=MCMCglmm(model.CBU.AOTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBU.CBL.2=MCMCglmm(model.CBU.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBL.AOTU.2=MCMCglmm(model.CBL.AOTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBL.CBU.2=MCMCglmm(model.CBL.CBU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.CBL.PB.2=MCMCglmm(model.CBL.PB, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.PB.POTU.2=MCMCglmm(model.PB.POTU, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
model.PB.CBL.2=MCMCglmm(model.PB.CBL, random=~phylo, family=c("gaussian"),
ginverse=list(phylo=inv.phylo$Ainv),
prior=prior, data=cx.data,
nitt=500000,burnin=10000,thin=500)
##diagnostics
## convergence
gelman.diag(mcmc.list(model.AOTU.CBU.1$Sol,model.AOTU.CBU.2$Sol))
gelman.diag(mcmc.list(model.AOTU.CBU.1$VCV,model.AOTU.CBU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.AOTU.CBU.1$VCV,model.AOTU.CBU.2$VCV))
plot(mcmc.list(model.AOTU.CBU.1$Sol,model.AOTU.CBU.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.AOTU.CBL.1$Sol,model.AOTU.CBL.2$Sol))
gelman.diag(mcmc.list(model.AOTU.CBL.1$VCV,model.AOTU.CBL.2$VCV))
#random and fixed okay
plot(mcmc.list(model.AOTU.CBL.1$VCV,model.AOTU.CBL.2$VCV))
plot(mcmc.list(model.AOTU.CBL.1$Sol,model.AOTU.CBL.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.POTU.PB.1$Sol,model.POTU.PB.2$Sol))
gelman.diag(mcmc.list(model.POTU.PB.1$VCV,model.POTU.PB.2$VCV))
#random and fixed okay
plot(mcmc.list(model.POTU.PB.1$VCV,model.POTU.PB.2$VCV))
plot(mcmc.list(model.POTU.PB.1$Sol,model.POTU.PB.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.CBU.AOTU.1$Sol,model.CBU.AOTU.2$Sol))
gelman.diag(mcmc.list(model.CBU.AOTU.1$VCV,model.CBU.AOTU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBU.AOTU.1$VCV,model.CBU.AOTU.2$VCV))
plot(mcmc.list(model.CBU.AOTU.1$Sol,model.CBU.AOTU.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.CBU.CBL.1$Sol,model.CBU.CBL.2$Sol))
gelman.diag(mcmc.list(model.CBU.CBL.1$VCV,model.CBU.CBL.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBU.CBL.1$VCV,model.CBU.CBL.2$VCV))
plot(mcmc.list(model.CBU.CBL.1$Sol,model.CBU.CBL.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.CBL.AOTU.1$Sol,model.CBL.AOTU.2$Sol))
gelman.diag(mcmc.list(model.CBL.AOTU.1$VCV,model.CBL.AOTU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.AOTU.1$VCV,model.CBL.AOTU.2$VCV))
plot(mcmc.list(model.CBL.AOTU.1$Sol,model.CBL.AOTU.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.CBL.CBU.1$Sol,model.CBL.CBU.2$Sol))
gelman.diag(mcmc.list(model.CBL.CBU.1$VCV,model.CBL.CBU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.CBU.1$VCV,model.CBL.CBU.2$VCV))
plot(mcmc.list(model.CBL.CBU.1$Sol,model.CBL.CBU.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.CBL.PB.1$Sol,model.CBL.PB.2$Sol))
gelman.diag(mcmc.list(model.CBL.PB.1$VCV,model.CBL.PB.2$VCV))
#random and fixed okay
plot(mcmc.list(model.CBL.PB.1$VCV,model.CBL.PB.2$VCV))
plot(mcmc.list(model.CBL.PB.1$Sol,model.CBL.PB.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.PB.POTU.1$Sol,model.PB.POTU.2$Sol))
gelman.diag(mcmc.list(model.PB.POTU.1$VCV,model.PB.POTU.2$VCV))
#random and fixed okay
plot(mcmc.list(model.PB.POTU.1$VCV,model.PB.POTU.2$VCV))
plot(mcmc.list(model.PB.POTU.1$Sol,model.PB.POTU.2$Sol))
#all good.
## convergence
gelman.diag(mcmc.list(model.PB.CBL.1$Sol,model.PB.CBL.2$Sol))
gelman.diag(mcmc.list(model.PB.CBL.1$VCV,model.PB.CBL.2$VCV))
#random and fixed okay
plot(mcmc.list(model.PB.CBL.1$VCV,model.PB.CBL.2$VCV))
plot(mcmc.list(model.PB.CBL.1$Sol,model.PB.CBL.2$Sol))
#all good.
sum.AOTU.CBU=summary(model.AOTU.CBU.1)
sum.AOTU.CBU.output=sum.AOTU.CBU$solutions
sum.AOTU.CBL=summary(model.AOTU.CBL.1)
sum.AOTU.CBL.output=sum.AOTU.CBL$solutions
sum.POTU.PB=summary(model.POTU.PB.1)
sum.POTU.PB.output=sum.POTU.PB$solutions
sum.CBU.AOTU=summary(model.CBU.AOTU.1)
sum.CBU.AOTU.output=sum.CBU.AOTU$solutions
sum.CBU.CBL=summary(model.CBU.CBL.1)
sum.CBU.CBL.output=sum.CBU.CBL$solutions
sum.CBL.AOTU=summary(model.CBL.AOTU.1)
sum.CBL.AOTU.output=sum.CBL.AOTU$solutions
sum.CBL.CBU=summary(model.CBL.CBU.1)
sum.CBL.CBU.output=sum.CBL.CBU$solutions
sum.CBL.PB=summary(model.CBL.PB.1)
sum.CBL.PB.output=sum.CBL.PB$solutions
sum.PB.POTU=summary(model.PB.POTU.1)
sum.PB.POTU.output=sum.PB.POTU$solutions
sum.PB.CBL=summary(model.PB.CBL.1)
sum.PB.CBL.output=sum.PB.CBL$solutions
#no pollen feeding effects
##### Rate analysis w multirateBM (phytools) (related to Figure 2 and S3) #####
setwd("C:/R-Analysis")
rate.data=read.table(file="Hel-CX-rates_vs2.txt", sep="\t", header=TRUE)
str(rate.data)
## import CX and rCBR data and reformat as named vectors
names(rate.data)[1] <- ""
rate.data <- data.frame(rate.data[,-1], row.names=rate.data[,1])
CX<-setNames(rate.data$log.CX,
rownames(rate.data))
rCBR<-setNames(rate.data$log.rCBR,
rownames(rate.data))
MB<-setNames(rate.data$log.MB,
rownames(rate.data))
AOTU<-setNames(rate.data$log.AOTU,
rownames(rate.data))
POTU<-setNames(rate.data$log.POTU,
rownames(rate.data))
## import tree
tt = read.nexus("Heliconiini.trees")
## Fit multirateBM
fit.HeliconiiniCX<-multirateBM(tt,CX,n.iter=3)
fit.HeliconiinirCBR<-multirateBM(tt,rCBR,n.iter=3)
fit.HeliconiiniMB<-multirateBM(tt,MB,n.iter=3)
fit.HeliconiiniAOTU<-multirateBM(tt,AOTU,n.iter=3)
fit.HeliconiiniPOTU<-multirateBM(tt,POTU,n.iter=3)
## CX
sig2 <- fit.HeliconiiniCX$sig2
fit.HeliconiiniCX$sig2 <- sig2
fit.HeliconiiniCX$tree <- tt
##rCBR
sig2 <- fit.HeliconiinirCBR$sig2
fit.HeliconiinirCBR$sig2 <- sig2
fit.HeliconiinirCBR$tree <- tt
##MB
sig2 <- fit.HeliconiiniMB$sig2
fit.HeliconiiniMB$sig2 <- sig2
fit.HeliconiiniMB$tree <- tt
##AOTU
sig2 <- fit.HeliconiiniAOTU$sig2
fit.HeliconiiniAOTU$sig2 <- sig2
fit.HeliconiiniAOTU$tree <- tt
##POTU
sig2 <- fit.HeliconiiniPOTU$sig2
fit.HeliconiiniPOTU$sig2 <- sig2
fit.HeliconiiniPOTU$tree <- tt
## calculate edge values function
ln.mean<-function(x){
if(x[1]==x[2]) return(x[1])
else {
a<-x[2]
b<-log(x[1])-log(x[2])
return(a/b*exp(b)-a/b)
}
}
##CX rates
CX.sig2<-apply(fit.HeliconiiniCX$tree$edge,1,function(e,fit.HeliconiiniCX)
ln.mean(fit.HeliconiiniCX[e]),fit=fit.HeliconiiniCX$sig2)
## rCBR rates
rCBR.sig2<-apply(fit.HeliconiinirCBR$tree$edge,1,function(e,fit.HeliconiinirCBR)
ln.mean(fit.HeliconiinirCBR[e]),fit=fit.HeliconiinirCBR$sig2)
## MB rates
MB.sig2<-apply(fit.HeliconiiniMB$tree$edge,1,function(e,fit.HeliconiiniMB)
ln.mean(fit.HeliconiiniMB[e]),fit=fit.HeliconiiniMB$sig2)
## AOTU rates
AOTU.sig2<-apply(fit.HeliconiiniAOTU$tree$edge,1,function(e,fit.HeliconiiniAOTU)
ln.mean(fit.HeliconiiniAOTU[e]),fit=fit.HeliconiiniAOTU$sig2)
## POTU rates
POTU.sig2<-apply(fit.HeliconiiniPOTU$tree$edge,1,function(e,fit.HeliconiiniPOTU)
ln.mean(fit.HeliconiiniPOTU[e]),fit=fit.HeliconiiniPOTU$sig2)
## regress CX sig2 against rCBR sig2 for branches and get residuals
CX.rCBR <- lm(CX.sig2~rCBR.sig2)
residuals.CX <- CX.rCBR$residuals
## regress MB sig2 against rCBR sig2 for branches and get residuals
MB.rCBR <- lm(MB.sig2~rCBR.sig2)
residuals.MB <- MB.rCBR$residuals
## regress AOTU sig2 against rCBR sig2 for branches and get residuals
AOTU.rCBR <- lm(AOTU.sig2~rCBR.sig2)
residuals.AOTU <- AOTU.rCBR$residuals
## regress POTU sig2 against rCBR sig2 for branches and get residuals
POTU.rCBR <- lm(POTU.sig2~rCBR.sig2)
residuals.POTU <- POTU.rCBR$residuals
## set residuals minimum to 0
residuals.CX <- residuals.CX - min(residuals.CX)
residuals.MB <- residuals.MB - min(residuals.MB)
residuals.AOTU <- residuals.AOTU - min(residuals.AOTU)
residuals.POTU <- residuals.POTU - min(residuals.POTU)
## plot the residuals along branches. CX according to MB rates.
plot<- TRUE
cols<-setNames(wes_palette("Zissou1", 1000, type = "continuous"),
1:1000)
min.sig2<-min(residuals.MB)
max.sig2<-max(residuals.MB)
edge.states<-vector()
for(i in 1:length(residuals.CX)){
edge.states[i]<-round((residuals.CX[i]-min.sig2)/
(max.sig2-min.sig2)*999)+1
}
tree<-paintBranches(fit.HeliconiiniCX$tree,edge=fit.HeliconiiniCX$tree$edge[1,2],
state=edge.states[1])
for(i in 2:length(edge.states))
tree<-paintBranches(tree,edge=tree$edge[i,2],
state=edge.states[i])
nticks<-10
ticks<-seq(min(residuals.CX),max(residuals.CX),
length.out=nticks)
object<-list(tree=tree,cols=cols,ticks=ticks)
class(object)<-"multirateBM_plot"
plot(object,outline=FALSE)
invisible(object)
## plot the residuals along branches. CX normal
plot<- TRUE
cols<-setNames(wes_palette("Zissou1", 1000, type = "continuous"),
1:1000)
min.sig2<-min(residuals.CX)
max.sig2<-max(residuals.CX)
edge.states<-vector()
for(i in 1:length(residuals.CX)){
edge.states[i]<-round((residuals.CX[i]-min.sig2)/
(max.sig2-min.sig2)*999)+1
}
tree<-paintBranches(fit.HeliconiiniCX$tree,edge=fit.HeliconiiniCX$tree$edge[1,2],
state=edge.states[1])
for(i in 2:length(edge.states))
tree<-paintBranches(tree,edge=tree$edge[i,2],
state=edge.states[i])
nticks<-10
ticks<-seq(min(residuals.CX),max(residuals.CX),
length.out=nticks)
object<-list(tree=tree,cols=cols,ticks=ticks)
class(object)<-"multirateBM_plot"
plot(object,outline=FALSE)
invisible(object)
## plot the residuals along branches. MB normal
plot<- TRUE
cols<-setNames(wes_palette("Zissou1", 1000, type = "continuous"),
1:1000)
min.sig2<-min(residuals.MB)
max.sig2<-max(residuals.MB)
edge.states<-vector()
for(i in 1:length(residuals.MB)){
edge.states[i]<-round((residuals.MB[i]-min.sig2)/
(max.sig2-min.sig2)*999)+1
}
tree<-paintBranches(fit.HeliconiiniMB$tree,edge=fit.HeliconiiniMB$tree$edge[1,2],
state=edge.states[1])
for(i in 2:length(edge.states))
tree<-paintBranches(tree,edge=tree$edge[i,2],
state=edge.states[i])
nticks<-10
ticks<-seq(min(residuals.MB),max(residuals.MB),
length.out=nticks)
object<-list(tree=tree,cols=cols,ticks=ticks)
class(object)<-"multirateBM_plot"
plot(object,outline=FALSE)
invisible(object)
## plot the residuals along branches. AOTU normal
plot<- TRUE
cols<-setNames(wes_palette("Zissou1", 1000, type = "continuous"),
1:1000)
min.sig2<-min(residuals.AOTU)
max.sig2<-max(residuals.AOTU)
edge.states<-vector()
for(i in 1:length(residuals.AOTU)){
edge.states[i]<-round((residuals.AOTU[i]-min.sig2)/
(max.sig2-min.sig2)*999)+1
}
tree<-paintBranches(fit.HeliconiiniAOTU$tree,edge=fit.HeliconiiniAOTU$tree$edge[1,2],
state=edge.states[1])
for(i in 2:length(edge.states))
tree<-paintBranches(tree,edge=tree$edge[i,2],
state=edge.states[i])
nticks<-10
ticks<-seq(min(residuals.AOTU),max(residuals.AOTU),
length.out=nticks)
object<-list(tree=tree,cols=cols,ticks=ticks)
class(object)<-"multirateBM_plot"
plot(object,outline=FALSE)
invisible(object)
## plot the residuals along branches. POTU normal
plot<- TRUE
cols<-setNames(wes_palette("Zissou1", 1000, type = "continuous"),
1:1000)
min.sig2<-min(residuals.POTU)
max.sig2<-max(residuals.POTU)
edge.states<-vector()
for(i in 1:length(residuals.POTU)){
edge.states[i]<-round((residuals.POTU[i]-min.sig2)/
(max.sig2-min.sig2)*999)+1
}
tree<-paintBranches(fit.HeliconiiniPOTU$tree,edge=fit.HeliconiiniPOTU$tree$edge[1,2],
state=edge.states[1])
for(i in 2:length(edge.states))
tree<-paintBranches(tree,edge=tree$edge[i,2],
state=edge.states[i])
nticks<-10
ticks<-seq(min(residuals.POTU),max(residuals.POTU),
length.out=nticks)
object<-list(tree=tree,cols=cols,ticks=ticks)
class(object)<-"multirateBM_plot"
plot(object,outline=FALSE)
invisible(object)
#plot rates
rate.data.mod=data.frame(CX.sig2, rCBR.sig2,residuals)
str(rate.data.mod)
rate.plot=ggplot(rate.data.mod,aes(x=rCBR.sig2, y=CX.sig2))+geom_point(alpha=0.4,size=4)
rate.plot
rate.plot2=rate.plot+ geom_smooth(method = "lm",size=3) + theme_minimal()
rate.plot2
ggplotly(rate.plot2)
summary(CX.rCBR)
##### ER neuron cell count (Figure 6) #####
setwd("C:/R-Analysis")
cellcount.data=read.table(file="cellcounts_.txt", sep="\t", header=TRUE)
str(cellcount.data)
hel.data=subset(cellcount.data, species == "hel")
iulia.data=subset(cellcount.data, species == "iulia")
count.lm.iul=lm(cell_count~no_sections, iulia.data)
count.lm.hel=lm(cell_count~no_sections, hel.data)
summary(count.lm.iul)
summary(count.lm.hel)
count.lm=lm(cell_count~no_sections, cellcount.data)
summary(count.lm)
#no of sections that were merged does not influence number of cells.
#cell numbers in GAD in cx
setwd("C:/R-Analysis")
count.data=read.table(file="Cell-count_CX.txt", sep="\t", header=TRUE)
str(count.data)
count.ER.lm=lm(number~spec+sex,data=count.data)
summary(count.ER.lm)
count.plot=ggplot(count.data, aes(x=spec,y=number,fill=spec))
count.plot
count.plot2=count.plot+theme_minimal()+geom_dotplot(binaxis='y', stackdir='center',stackratio=1.5, dotsize=1.2)
count.plot2
count.plot3=count.plot2+scale_fill_manual(values = c("#1293b9ff", "#cb8781ff"))
count.plot3
##### bulb quantification (Figure 7) #####
setwd("C:/R-Analysis")
bulb.data=read.table(file="Bulb_quant.txt", sep="\t", header=TRUE)
str(bulb.data)
bulb.pers=lm(log10_bulb.vol~bulb.MF.TH+spec+sex,data=bulb.data)
MG.pers=lm(MG.AVG~MG.MF.TH+spec+sex,data=bulb.data)
summary(bulb.pers)
summary(MG.pers)
#person doing segmentations does not harbour specific information.
#are there species and sex differences? we dont expect these two factors to interact in a two-species comparison.
bulb.spec=lm(log10_bulb.vol~spec+sex,data=bulb.data)
summary(bulb.spec)
#no species nor sex differences.
#is total MG no predicted by avg mg size, or bulb volume?
total.MG.no=lm(No_MG.bulb~log10_bulb.vol+MG.AVG+spec+sex,data=bulb.data)
summary(total.MG.no)
#bulb volume and avg MG are predictors of total MG no, which makes sense, as this a pure mathematical value and derived from them.
#does ER neuron number predict bulb vol, or total MG no, or avg MG size (which is a reasonable assumption as neuron number might predict these metrics)?
bulb.ER=lm(log10_bulb.vol~ER.no+spec+sex,data=bulb.data)
summary(bulb.ER)
MGno.ER=lm(No_MG.bulb~ER.no+spec+sex,data=bulb.data)
summary(MGno.ER)
MGavg.ER=lm(MG.AVG~ER.no+spec+sex,data=bulb.data)
summary(MGavg.ER)
#none of them do. So, it seems like ER neuron number is largely independent in this case, from these bulb related metrics.
#plotting
bulb.plot=ggplot(bulb.data, aes(x=spec,y=log10_bulb.vol,fill=spec))
bulb.plot
bulb.plot2=bulb.plot+theme_minimal()+geom_dotplot(binaxis='y', stackdir='center',stackratio=1.5, dotsize=1.2)
bulb.plot2
bulb.plot3=bulb.plot2+scale_fill_manual(values = c("#1293b9ff", "#cb8781ff"))
bulb.plot3
MGno.plot=ggplot(bulb.data, aes(x=spec,y=No_MG.bulb,fill=spec))
MGno.plot
MGno.plot2=MGno.plot+theme_minimal()+geom_dotplot(binaxis='y', stackdir='center',stackratio=1.5, dotsize=1.2)
MGno.plot2
MGno.plot3=MGno.plot2+scale_fill_manual(values = c("#1293b9ff", "#cb8781ff"))
MGno.plot3
#to plot MG sizes
setwd("C:/R-Analysis")
MG_avg.data=read.table(file="MG.avg.txt", sep="\t", header=TRUE)
str(MG_avg.data)
MGsize.plot=ggplot(MG_avg.data, aes(x=spec,y=MG.size,fill=spec))
MGsize.plot
MGsize.plot2=MGsize.plot+theme_minimal()+geom_dotplot(binaxis='y', stackdir='center',stackratio=1.5, dotsize=1.2)
MGsize.plot2
MGsize.plot3=MGsize.plot2+scale_fill_manual(values = c("#1293b9ff", "#cb8781ff"))
MGsize.plot3
bulb.mg.plot=ggarrange(MGno.plot3,MGsize.plot3,bulb.plot3,
ncol = 3, nrow = 1,common.legend = TRUE, legend="bottom")
bulb.mg.plot

Heliconiini-CX_SCRIPT.R, under CC-BY-4.0 · at the source

Overview

Authors: Max S Farnworth¹, Yi Peng Toh², Theodora Loupasaki¹, Elizabeth A Hodge¹, Basil el Jundi^3,4, Stephen H Montgomery¹

ORCID iDs: Max S Farnworth, Yi Peng Toh, Theodora Loupasaki, Elizabeth A Hodge, Basil el Jundi, Stephen H Montgomery

Institutions: University of Bristol (United Kingdom); Ludwig-Maximilians-Universität München (Germany); Carl von Ossietzky Universität Oldenburg (Germany); Norwegian University of Science and Technology (Norway)

Journal: eLife, volume 14, article RP107589

Dates: published online 1 June 2026

Type: Research article · Language: English

License: CC BY

Identifiers: DOI 10.7554/elife.107589 · PMID 42223014 · PMCID PMC13225844 · OpenAlex W4414262609

Open access: gold, a free copy (OpenAlex)

Status: code verified

Categories: other (organism)

Methods: Statistics

Keywords: Other

MeSH: Biological Evolution*, Butterflies*, Mushroom Bodies*, Animals, Neuropeptides (* major topic)

Topic: Plant and animal studies (Ecology, Evolution, Behavior and Systematics, Agricultural and Biological Sciences), according to OpenAlex

Funding: European Research Council (Starter Grant 851040, to Richard M Merrill, 851040, Starter Grant 758508, 758508); Natural Environment Research Council (Independent Research Fellowship NE/N014936/1); Human Frontier Science Program (Project Grant RGP0029/2022); Medical Research Council (DTP studentship MR/W006308/1); University of Bristol (Student Scholarship); Deutsche Forschungsgemeinschaft (Walter-Benjamin Fellowship FA 1818/1-1)

Citations: cited by 1 paper (Europe PMC); 110 references in the paper

Research resources: Anti-FASII (mouse; monoclonal) RRID:AB_10805878, Anti-mouse-Cyanine-2 (goat; polyclonal) RRID:AB_2307343, Anti-rabbit-Cyanine-3 (goat; polyclonal) RRID:AB_2338006, Anti-rabbit-Cyanine-2 (goat; polyclonal) RRID:AB_2338021, Anti-mouse-Cyanine-3 (goat; polyclonal) RRID:AB_2338680, Anti-HRP (rabbit; polyclonal) RRID:AB_261181, Anti-TH (rabbit; polyclonal) RRID:AB_390204, Anti-GAD (rabbit; polyclonal) RRID:AB_477019, Anti-SYNORF1 (mouse; monoclonal) RRID:AB_528479, Anti-5-HT (rabbit; polyclonal) RRID:AB_572263, Anti-acTUB (mouse; monoclonal) RRID:AB_609894, R Studio 2023.06.0 RRID:SCR_000432, R 4.3.1 RRID:SCR_001905, Fiji 1.54 RRID:SCR_003070, Amira 2020.2/5.4.3/2021.1 RRID:SCR_007353

Abstract

Neural circuits evolved to produce variable cognitive processes through adaptive mechanisms operating within a background of developmental and functional constraints. Understanding how this conflict is resolved requires a comparative framework encapsulating clear behavioural variation. We leverage Heliconiini butterflies to examine how selection shaped the evolution of the central complex and mushroom bodies, two insect integration centres involved in navigation. The evolution of systematic spatial foraging in Heliconius has led to changes in brain morphology and learning and memory profiles over a short evolutionary timescale. Here, we show that in contrast to massively expanded mushroom bodies, the central complex is strongly conserved in size and general architecture. However, we identify divergences in the expression of a neuropeptide, Allatostatin A, in the noduli, and in the numbers of GABA-ergic ring neurons and their branching in the fan-shaped body, which are essential members of the anterior compass pathway. These differences are rare examples of divergence inside the central complex network matching expectations of where evolutionary adaptability might occur. We conclude that due to the contrasting volumetric conservation of the central complex, and the massive differences in the mushroom bodies, their circuit logics must determine distinct responses to selection associated with divergent foraging behaviours.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 3 matches between paragraphs and lines of code.

Zenodo 15304965

License: CC-BY-4.0

State: the link answers, verified on 27 September 2026

Evidence: files inventoried

Languages: R (1)

Size: 10 files, 1 script

Software Heritage: not checked

Found in: the text, “Statistical tests of volumetric dataset”

Not found: README, license file, CITATION.cff, environment file, tests, continuous integration, documentation

Tools: emmeans (1 file), ggpubr (1 file), Plotly (1 file), tidyverse (1 file)

Availability: 1 check, the latest on 27 September 2026: the link answers (HTTP 200)

27 September 2026: the link answers (HTTP 200)

1 file

Heliconiini-CX_SCRIPT.R, R, 2,308 lines, 3 matches

At the source: doi.org/10.5281/zenodo.15304965

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
1 script, each with its path and the digest of its content;
3 matches between paragraphs of the paper and lines of the code (method lexical-v1);
neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

Data, R script and readme files associated with this publication is publicly available at https://doi.org/10.5281/zenodo.15304965.

The following dataset was generated:

Farnworth MS. 2025. Dataset associated with publication of Farnworth et al "Distinct evolutionary trajectories of two integration centres, the central complex and mushroom bodies, across Heliconiini butterflies". Zenodo.

The following previously published dataset was used:

Couto A, Young F, Atzeni D, Marty S, Melo-Flórez L, Hebberecht L, Monllor M, Neal C, Cicconardi F, McMillan W, Montgomery SH. 2023. Data From: Rapid expansion and visual specialisation of learning and memory centers in the brains of Heliconiini butterflies. Dryad Digital Repository.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 27 September 2026: the first record

Recorded: type, language, journal, volume, pages, dates, 6 authors, 1 keyword, 5 MeSH terms, 6 funders, 108 references, 15 RRIDs.

Cite

This paper

Farnworth, M. S., Toh, Y. P., Loupasaki, T., Hodge, E. A., el Jundi, B., & Montgomery, S. H. (2026). Distinct evolutionary trajectories of two integration centres, the central complex and mushroom bodies, across Heliconiini butterflies. eLife, 14, RP107589. https://doi.org/10.7554/elife.107589

BibTeX

@article{farnworth2026distinct, author = {Farnworth, Max S and Toh, Yi Peng and Loupasaki, Theodora and Hodge, Elizabeth A and el Jundi, Basil and Montgomery, Stephen H}, title = {{Distinct evolutionary trajectories of two integration centres, the central complex and mushroom bodies, across Heliconiini butterflies}}, journal = {eLife}, year = {2026}, month = jun, volume = {14}, pages = {RP107589}, publisher = {eLife Sciences Publications, Ltd}, issn = {2050-084X}, doi = {10.7554/elife.107589}, url = {https://doi.org/10.7554/elife.107589}, pmid = {42223014}, pmcid = {PMC13225844} }

RIS

TY - JOUR AU - Farnworth, Max S AU - Toh, Yi Peng AU - Loupasaki, Theodora AU - Hodge, Elizabeth A AU - el Jundi, Basil AU - Montgomery, Stephen H TI - Distinct evolutionary trajectories of two integration centres, the central complex and mushroom bodies, across Heliconiini butterflies T2 - eLife J2 - eLife PY - 2026 DA - 2026/06/01 VL - 14 SP - RP107589 SN - 2050-084X PB - eLife Sciences Publications, Ltd DO - 10.7554/elife.107589 UR - https://doi.org/10.7554/elife.107589 LA - en ER -

CSL-JSON

{ "id": "10.7554/elife.107589", "type": "article-journal", "title": "Distinct evolutionary trajectories of two integration centres, the central complex and mushroom bodies, across Heliconiini butterflies", "container-title": "eLife", "author": [ { "family": "Farnworth", "given": "Max S" }, { "family": "Toh", "given": "Yi Peng" }, { "family": "Loupasaki", "given": "Theodora" }, { "family": "Hodge", "given": "Elizabeth A" }, { "family": "el Jundi", "given": "Basil" }, { "family": "Montgomery", "given": "Stephen H" } ], "container-title-short": "eLife", "volume": "14", "page": "RP107589", "DOI": "10.7554/elife.107589", "PMID": "42223014", "PMCID": "PMC13225844", "ISSN": "2050-084X", "publisher": "eLife Sciences Publications, Ltd", "URL": "https://doi.org/10.7554/elife.107589", "language": "en", "issued": { "date-parts": [ [ 2026, 6, 1 ] ] } }

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.