# Variation in the strength of allometry drives rates of evolution in primate brain shape # Package library(RRphylo) library(ape) library(ggplot2) library(ggtree) library(geiger) library(EBImage) library(LambertW) ### load workspace load(file="..Path/Workspace.Rdata") ## Analyses ### PCA and Shape analysis. To visualize the associated shape chages open the html file "Shape_changes.html" ### plot(match.scores[,c(1,2)], ylab="PC2",xlab="PC1",cex=2,mgp=c(1.5,0.5,0),font.lab=2,pch=21,col="black",bg="gray", asp=1, xlim=range(scores[,1])*1.3,ylim=range(scores[,2])*1.3) for(i in 1:length(unique(newf))){ subset(match.scores,as.character(newf)==unique(as.character(newf))[i])->xyscs points(xyscs[,1:2], cex=2,pch=21,col="black",bg=unique(as.numeric(newf))[i]) Plot_ConvexHull(xcoord = xyscs[,1], ycoord = xyscs[,2], lcolor = unique(as.numeric(newf))[i]) } legend("topright",legend = unique(as.character(newf)), fill=unique(as.numeric(newf)), bty="n", cex=0.8, x.intersp = 0.3) ## Multivariate Regressions Shape vs Size ccaprim<-mycca1(log(ss),matcha) plot(log(ss),ccaprim$yscores, ylab="Procrustes coordinates",xlab="ln Centroid Size", cex=2,mgp=c(1.5,0.5,0),font.lab=2,pch=21,col="black",bg="gray", asp=1, axes=F, ylim = range(ccaprim$yscores)*1.5) axis(1, pos=0, c(5,5.5,6,6.5,7,7.5)) axis(2, pos=6.25, c(-1,-0.8,-0.4,0, 0.4,0.8,1)) for(i in 1:length(unique(newf))){ subset(data.frame(log(ss),ccaprim$yscores),as.character(newf)==unique(as.character(newf))[i])->xyscs points(xyscs[,1:2], cex=2,pch=21,col="black",bg=unique(as.numeric(newf))[i]) Plot_ConvexHull(xcoord = xyscs[,1], ycoord = xyscs[,2], lcolor = unique(as.numeric(newf))[i]) } legend("topleft",legend = unique(as.character(newf)), fill=unique(as.numeric(newf)), bty="n", cex=0.75, x.intersp = 0.2 ) ### RRphylo analyses ## Compute evolutionary rates prim.rates.p<-RRphylo(tree=tro,y=scores) # search shifts shifts.a.p<-search.shift(prim.rates.p,status.type = 'clade', foldername = "/Volumes/TOSHIBA EXT/Primates R codes/Meshes") ## visualize rates of focal clades dev.off() plotRates(prim.rates.p,node = 272,export.tiff = FALSE) plotRates(prim.rates.p,node = 231,export.tiff = c(FALSE)) plotRates(prim.rates.p,node = 213,export.tiff = c(FALSE)) plotRates(prim.rates.p,node = 170,export.tiff = c(FALSE)) ## test signficance in focal clades shifts.a.sing.p<-search.shift(prim.rates.p,status.type = 'clade',node=c(272,231,213,170), foldername = "/Volumes/TOSHIBA EXT/Primates R codes/Meshes") ## test statistics shifts.a.sing.p$single.clades ### Visualize shifts on the phylogenetic tree ggtree(tro, layout="fan")->p p + #geom_tiplab(size=3) geom_hilight(node=272, fill="red", alpha=0.2)+ geom_hilight(node=230, fill="red", alpha=0.2)+ geom_hilight(node=212, fill="navyblue", alpha=0.2)+ geom_hilight(node=170, fill="navyblue", alpha=0.2)+ geom_strip(taxa1 = "Pongo_abelii",taxa2 = "Hylobates_klossii",barsize=5, color="#6379ac",offset = 2)+ geom_strip(taxa1 = "Colobus_polykomos",taxa2 = "Semnopithecus_schistaceus",barsize=5, color="#4cc185",offset = 2)+ geom_strip(taxa1 = "Allenopithecus_nigroviridis",taxa2 = "Cercopithecus_mitis",barsize=5, color="#666363",offset = 2)+ geom_strip(taxa1 = "Macaca_fascicularis",taxa2 = "Papio_hamadryas",barsize=5, color="#f26065",offset = 2)+ geom_strip(taxa1 = "Nycticebus_coucang",taxa2 = "Euoticus_elegantulus",barsize=5, color="#fc6600",offset = 2)+ geom_strip(taxa1 = "Propithecus_verreauxi",taxa2 = "Varecia_variegata_variegata",barsize=5, color="#fcf415",offset = 2)+ geom_strip(taxa1 = "Cacajao_calvus",taxa2 = "Mico_argentatus",barsize=5, color="#afe4fa",offset = 2)+ geom_nodepoint(aes(color=cols),size=5)+ scale_color_gradient(low = "cyan1",high = "magenta2")+ geom_nodepoint(shape=21,size=5)+ geom_cladelabel(node=271, label="Strepsirrhini",hjust = "center", angle=-57,offset.text=7.5, barsize=0,color='black',offset = 2.5,fontsize = 4) + geom_cladelabel(node=212, label="Hominoidea", angle=285,hjust = "center",offset.text=7.5,barsize=0,color = 'black',offset = 2.5,fontsize = 4) + geom_cladelabel(node=170, label="Papionini", angle=65,hjust = "center",offset.text=7.5,barsize=0,color='black',offset = 2.5,fontsize = 4) + geom_cladelabel(node=230, label="Platyrrhini", angle=-325,hjust = "center",offset.text=7.5,barsize=0,color='black',offset = 2.5,fontsize = 4) + geom_cladelabel(node=189, label="Colobinae", angle=-25,hjust = "center",offset.text=7.5,barsize=0,color='black',offset = 2.5,fontsize = 4) + geom_cladelabel(node=157, label="Cercopithecini", angle=25,hjust = "center",offset.text=7.5,barsize=0,color='black',offset = 2.5,fontsize = 4) #scale_color_gradientn(colours = c("red","blue")) theme(legend.position = "none") ### visualize angles and trajectories difference dev.off() datu<-ggplot(angu.1, aes(x=angi.1, y=veci.1)) pu <- datu + theme_bw() pu+ geom_segment(aes(y=0, xend = angi.1, yend = veci.1, color = V1), size = 1.5, alpha = 0.7) + geom_vline(xintercept = seq(0, 360-1, by = 45), colour = "grey90", size = 0.2) + coord_polar(theta='x',start = 11,direction=1,clip='on')+ labs(x = 'Angles', y = 'Vector Size') + scale_x_continuous(limits = c(0, 360), expand = c(0, 0), breaks = seq(0, 360-1, by = 45)) ##