#qindex Rank #Uses the gplots package library(gplots) qset=read.table("~/Dropbox/The boss/Letocoris_ms/k2_structure_Qoutput",header=T) attach(qset) #orders the data by q value of Atalaya (-, in descending order) qset_ord=qset[order(-Qatalaya),] #adds a colunm with the rank of the 89 ordered samples qset_ord$qrank=c(1:89) #adds colunm to plot the uppperbar of the credible interval qset_ord$uperr=qset_ord$upAt-qset_ord$Qatalaya #adds colunm to plot the lowerbar of the credible interval qset_ord$lowerr=qset_ord$Qatalaya-qset_ord$lowAt #defines subset of ordered samples in C. halicacabum hali=subset(qset_ord, pop_ID==1 |pop_ID==2 |pop_ID==3 ) #defines the subset of ordered samples of the sympatric populations of Atalaya symp=subset(qset_ord, pop_ID==9 |pop_ID==10 ) #defines the subset of ordered samples of allopatric Atalaya AtaHost=subset(qset_ord, pop_ID==4 |pop_ID==5 |pop_ID==6 |pop_ID==7 |pop_ID==8) #defines the subset of ordered samples of the East coast (L.t. tagalicus) EastHost=subset(qset_ord, pop_ID==11 |pop_ID==12 |pop_ID==13 |pop_ID==14 |pop_ID==15|pop_ID==16|pop_ID==17) #defines frame for plotting par(oma=c(3, 7, 3, 3)) #plots C.halicacabum populations plotCI(hali$qrank,hali$Qatalaya,pch=22, cex=1.5, pt.bg="dark grey", uiw=hali$uperr,liw=hali$lowerr,gap=0, xaxt="n", xlim=c(0,89),ylim=c(0,1), xlab="",ps= 20, ylab="",bty="l",lty=1, sfrac=0.0001, barcol="dark grey", xaxt="n", slty=2, scol="grey") # adds sympatric Atalaya population to plot_1 plotCI(symp$qrank,symp$Qatalaya,pch=22,cex=1.5, pt.bg="white", uiw=symp$uperr,liw=symp$lowerr,gap=0, xaxt="n", xlab="",ylab="",bty="l",lty=1, sfrac=0.0001, barcol="dark grey",xaxt="n",slty=2,scol="grey", add=T) # adds sympatric Atalaya population to plot_2 (adds little symbol) plotCI(symp$qrank,symp$Qatalaya,pch=8,cex=.5, pt.bg="white", uiw=symp$uperr,liw=symp$lowerr,gap=0, xaxt="n", xlab="",ylab="",bty="l",lty=1, sfrac=0.0001, barcol="dark grey",xaxt="n",slty=2,scol="grey", add=T) # adds allopatric Atalaya populations to plot plotCI(AtaHost$qrank,AtaHost$Qatalaya,pch=22,cex=1.5,pt.bg="white", uiw=AtaHost$uperr,liw=AtaHost$lowerr,gap=0, xaxt="n", xlab="",ylab="",bty="l",lty=1, sfrac=0.0001, barcol="dark grey",xaxt="n",slty=2,scol="grey",add=T) # adds East populations to plot plotCI(EastHost$qrank,EastHost$Qatalaya,pch=21,cex=1.5, pt.bg="black", uiw=EastHost$uperr,liw=EastHost$lowerr,gap=0, xaxt="n", xlab="",ylab="",bty="l",lty=1, sfrac=0.0001, barcol="dark grey",xaxt="n",slty=2,scol="grey",add=T) mtext(expression(bold("Admixture index (q)")),side=2,line=3,cex=1.7,outer="T") #LInear Mixed model library(lme4) m=read.table("~/Dropbox/The boss/Letocoris_ms/Reducedhost_morphology_ind_measurements_R",header=T) attach(m) names(m) foo= lmer(beak~host+gender+(1|pop)) summary(foo) mcmcfoo=mcmcsamp(foo, n=10000) HPDinterval(mcmcfoo@fixef, 0.95) densityplot(mcmcfoo@fixef["(Intercept)",1:10000]) densityplot(mcmcfoo@fixef["(genderM",1:10000]) NumbNeg = subset(mcmcfoo@fixef["hostCardiospermum grandiflorum",1:10000], mcmcfoo@fixef["hostAtalaya hemiglauca",1:10000] < 0) 2*(length(NumbNeg)/10000) # Morpholgy analysis #1 library(ggplot2) library(gplots) library(gridExtra) allbeaks=read.table("~/Dropbox/The boss/Letocoris_ms/morphdata_4hosts_Reducedata_R",header=T) attach(allbeaks) p1=qplot(host,Mean_beak,data=allbeaks,geom="boxplot",xlab="", ylab="Beak length (mm)") + facet_grid(.~sex,scales="free", space="free") +scale_y_continuous(breaks=seq(4.75, 8, 0.50)) +theme(strip.text.x = element_text(size=14,face="bold"),axis.title.y = element_text(face="bold", colour="black", size=14, vjust=0.3), axis.ticks = element_blank(), axis.text.x = element_blank(),axis.text.y = element_text(colour="black", size=11), panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank(),panel.background=element_blank(), panel.grid.major.x=element_blank(), panel.grid.minor.x=element_blank(), panel.border = element_rect(linetype = "solid", colour = "black",size=1, fill=NA), axis.line = element_line(colour="black"), axis.ticks.y=element_line(size=2)) +scale_x_discrete(limits=c("A. tomentosus", "C. grandiflorum", "C. halicacabum", "A. hemiglauca")) + geom_vline(aes(xintercept=2.5), colour="black", linetype="dashed") p2=qplot(host,Mean_ratio,data=allbeaks,geom="boxplot",xlab="", ylab="Beak/body ratio") + facet_grid(.~sex,scales="free", space="free") +theme(strip.text.x = element_blank(),strip.background= element_blank(),axis.text.x = element_text(face="italic", colour="black", size=12),axis.title.y = element_text(face="bold", colour="black", size=14,vjust=0.3), axis.text.y = element_text(colour="black", size=11), panel.grid.minor.y=element_blank(), panel.grid.major.y=element_blank(), panel.background=element_blank(), panel.grid.major.x=element_blank(), panel.grid.minor.x=element_blank(),panel.border = element_rect(linetype = "solid", colour = "black",size=1, fill=NA),axis.line = element_line(colour="black"), axis.ticks.y=element_line(colour="black",size=1)) +scale_x_discrete(limits=c("A. tomentosus", "C. grandiflorum", "C. halicacabum", "A. hemiglauca")) + geom_vline(aes(xintercept=2.5), colour="black", linetype="dashed") grid.arrange(p1,p2,nrow=2) #2 library(ggplot2) library(gridExtra) morpho=read.table("~/Dropbox/The boss/Letocoris_ms/morphology_AlloSymp_data.txt", header=T) attach(morpho) fem=subset(morpho,sex=="Females") p1=qplot(host,beak,data=fem,geom="boxplot",xlab="", ylab="Beak length (mm)") +theme(axis.title.y = element_text(face="bold", colour="black", size=15, vjust=0.4),axis.line = element_line(colour="black"), panel.background=element_blank(), panel.grid.major.x=element_blank(), panel.grid.major.y=element_blank(), panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text.x = element_blank(), axis.text.y = element_text(colour="black", size=13), axis.title.x = element_blank(),panel.border = element_rect(linetype = "solid", colour = "black",size=1, fill=NA)) + geom_vline(aes(xintercept=3.5), colour="black", linetype="dashed") +scale_x_discrete(limits=c("Ata_a","Ata_s","hali","east")) +scale_y_continuous(breaks=seq(4.5,8.7, 1)) p2=qplot(host,whole,data=fem,geom="boxplot",xlab="", ylab="Body length (mm)") +theme(axis.title.y = element_text(face="bold", colour="black", size=15, vjust=0.4),axis.line = element_line(colour="black"), panel.background=element_blank(), panel.grid.major.x=element_blank(), panel.grid.major.y=element_blank(), panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text.x = element_blank(), axis.text.y = element_text(colour="black", size=13), axis.title.x = element_blank(),panel.border = element_rect(linetype = "solid", colour = "black",size=1, fill=NA)) + geom_vline(aes(xintercept=3.5), colour="black", linetype="dashed") +scale_x_discrete(limits=c("Ata_a","Ata_s","hali","east")) p3=qplot(host,beak/whole,data=fem,geom="boxplot",xlab="", ylab="beak/body ratio") +theme(axis.title.y = element_text(face="bold", colour="black", size=15, vjust=0.4),axis.line = element_line(colour="black"), panel.background=element_blank(), panel.grid.major.x=element_blank(), panel.grid.major.y=element_blank(), panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text.x = element_text(colour="black",size=12,face="italic"), axis.text.y = element_text(colour="black", size=13), axis.title.x = element_blank(),panel.border = element_rect(linetype = "solid", colour = "black",size=1, fill=NA)) +geom_vline(aes(xintercept=3.5), colour="black", linetype="dashed")+scale_x_discrete(limits=c("Ata_a","Ata_s","hali","east"),breaks=c("Ata_a","Ata_s","hali","east"), labels=c("A.hemiglauca", "A.hemiglauca", "C.halicacabum","Alectryon/C.grandiflorum")) +scale_y_continuous(breaks=seq(0.4,0.8, 0.1)) grid.arrange(p1,p2,p3, ncol=1) #morphology _Vs_Admixtureindex library(ggplot2) library(gridExtra) library(gplots) tru=read.table("~/Dropbox/The boss/Letocoris_ms/morph_Qindex_SympatrydataWhole", header=T) attach(tru) tru2=ggplot(tru, aes(x=q, y=ratio, shape=group)) + geom_point(cex=3.5) +coord_cartesian(ylim=c(0.4, 0.8)) +scale_y_continuous(breaks=seq(0.5,0.8, 0.1)) +theme_bw() +scale_shape_manual(values=c(1,19)) +geom_smooth(method=lm,se=F, linetype="dashed", colour="dark grey") +theme(legend.position="bottom", legend.title = element_blank(),legend.text = element_text(colour="black", size = 14, face = "italic"), axis.line = element_line(colour="black"), panel.border=element_blank(), panel.grid.major.x=element_blank(), panel.grid.major.y=element_blank(), panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text.x = element_text(colour="black", size=13), axis.text.y = element_text(colour="black", size=13), axis.title.y = element_text(face="bold", colour="black", size=14,vjust=0.2), axis.title.x = element_text(face="bold", colour="black", size=14,vjust=0.3), legend.key=element_blank()) +labs(x="Admixture index (q)", y= "Beak/body ratio") tru1=ggplot(tru, aes(x=q, y=beak, shape=group)) + geom_point(cex=3.5) +coord_cartesian(ylim=c(4, 8.5)) +scale_y_continuous(breaks=seq(4.5, 8.75, 1)) +theme_bw() +scale_shape_manual(values=c(1,19)) +geom_smooth(method=lm,se=F, linetype="dashed", colour="dark grey") +theme(legend.position="none", axis.line = element_line(colour="black"), panel.border=element_blank(), panel.grid.major.x=element_blank(), panel.grid.major.y=element_blank(), panel.grid.minor.x=element_blank(), panel.grid.minor.y=element_blank(), axis.text.x = element_blank(), axis.text.y = element_text(colour="black", size=13), axis.title.y = element_text(face="bold", colour="black", size=14,vjust=0.2), axis.title.x = element_blank()) +labs(x="Admixture index (q)", y= "Beak length (mm)") grid.arrange(tru1,tru2,ncol=1)