########################################## ############ Step 1. Get set up to work ######## ########################################## library(ENmisc) library(betareg) ########################################## ############ Step 2. Import sister.summary.table ######## ########################################## sis <-read.csv ("~/Dropbox/Geography of Selfing Review/data sisters/sister.geo.area.csv") sis <-read.csv ("sister.summary.table.csv") names(sis) sis$mate.diff <- factor(sis$mate.diff) ########################################## ############ Step 3A-D. subset so that species only appear once in the data set ######## ########################################## #3A stop and think::: make sure data.frame is sorted by genus and then by posterior probability high to low!!! #3B bind genus-spA, gen.spB sis$gen.A=paste(sis$genus,sis$sp.A,sep="_") sis$gen.B=paste(sis$genus,sis$sp.B,sep="_") #3C sis.subA=sis[!duplicated(sis$gen.A),] #3D get list of gen.A and gen.B and check for duplicates (this is boneheaded) tmp=append(sis.subA$gen.A,sis.subA$gen.B) dup.row=duplicated(tmp) dup.row=dup.row[220:438] sis.subA$dup.row=dup.row sis.unique=subset(sis.subA, dup.row=="FALSE") nrow(sis.unique) sis=sis.unique sis$mate.diff <- factor(sis$mate.diff) table(sis$mate.diff, sis$genus) nrow(sis) table(sis$mate.diff) ########################################## ############ Step 4. Add variables for later ######## ########################################## sis$ln.dist=log(sis$dist) #log transformed branch length ########################################## ############ Step 5. Does divergence time vary by sister pair mating system? ######## ########################################## m1= aov(ln.dist ~ mate.diff, weights=pp, data=sis) summary(m1) TukeyHSD(m1) sis$fitted = fitted.values(m1) s.s=subset(sis,(mate.diff=="s_s")) s.o=subset(sis,(mate.diff=="s_o")) o.o=subset(sis,(mate.diff=="o_o")) exp(mean(o.o$fitted)) exp(mean(s.o$fitted)) exp(mean(s.s$fitted)) exp(mean(sis$fitted)) tiff("AJB.Fig1.tiff", height = 3.5, width = 3.5, units = 'in', res=800) par(mar=c(2.5,4.3,1,1), mfrow=c(1,1),cex.lab=1, oma=c(0,0,0,0)) wtd.boxplot(sis$ln.dist~sis$mate.diff,weights=sis$pp,ylab=expression(paste("Divergence time (mya)")), main="",outcol="white",names=c("","",""), ylim=c(-2,4.2),boxwex=0.5,border = c("black","black","black"),yaxt='n', col=c("gray30","red","pink"),cex.lab=1.3) mtext("A", side=1, line=-12.8, at=1,cex=0.7) mtext("B,C", side=1, line=-13.75, at=2,cex=0.7) mtext("A,C", side=1, line=-12.4, at=3,cex=0.7) axis(side=1, at=c(1,2,3),labels=c("o-o","s-o","s-s")) axis(side=2, at=c(-1.05,0,1.1,2.1,3,3.91), labels=c("0.35","1","3","8","20","50"), cex.axis=0.8) dev.off() ########################################## ############ Step 6. Does co-occurrence vary by sister pair mating system? ######## ########################################## #transform co-occurence response variable to eliminate 0 and 1 ala Smithson and Verkuilen 2006 # (y · (n − 1) + 0.5) /n where n is the sample size sis$trans.overlap1 = (sis$overlap.min.res1 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(sis$trans.overlap1 ~ mate.diff , weights=pp,data=sis, type="BC") summary(m1) sis$trans.overlap1 = (sis$overlap.min.res0.5 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(sis$trans.overlap1 ~ mate.diff , weights=pp,data=sis, type="BC") summary(m1) sis$trans.overlap1 = (sis$overlap.min.res0.1 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(sis$trans.overlap1 ~ mate.diff , weights=pp, data=sis, type="BC") summary(m1) sis$trans.overlap1 = (sis$overlap.min.res0.05 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(sis$trans.overlap1 ~ mate.diff, weights=pp , data=sis, type="BC") summary(m1) sis$fitted = fitted.values(m1) s.s=subset(sis,(mate.diff=="s_s")) s.o=subset(sis,(mate.diff=="s_o")) o.o=subset(sis,(mate.diff=="o_o")) mean(s.s$fitted) mean(s.o$fitted) mean(o.o$fitted) ################################################### ### Step 7. do single genus deletions impact step 5 significance? ################################################### gens=as.character(unique(sort(sis$genus))) #genus list #make lists for p values P0.1 = c() P0.05 = c() for (i in 1:length(gens)) { sis.sub=subset(sis, !(genus==gens[i])) sis.sub$genus=factor(sis.sub$genus) sis.sub$trans.overlap1 = (sis.sub$overlap.min.res0.1 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(trans.overlap1 ~ mate.diff ,weights=pp , data=sis.sub, type="BC") sum=summary(m1) my.p0.1 = round(sum$coefficients$mean[3,4],digits=4) #p value s-o sis.sub$trans.overlap1 = (sis.sub$overlap.min.res0.05 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(trans.overlap1 ~ mate.diff , weights=pp ,data=sis.sub, type="BC") sum=summary(m1) my.p0.05 = round(sum$coefficients$mean[3,4],digits=4) #p value s-o P0.1=c(P0.1,my.p0.1) P0.05=c(P0.05,my.p0.05) } mytable = cbind(P0.1,P0.05) colnames(mytable) = c("ss.p0.1", "ss.p0.05") #note 0.1 and 0.05 refers to the raster cell size in decimal degrees rownames(mytable) = gens mytable ################################################### ### Step 8. Does co-occurrence vary by sister pair age? ################################################### sis$trans.overlap1 = (sis$overlap.min.res1 * (nrow(sis) - 1) + 0.5) / nrow(sis) gy_logit <- betareg(trans.overlap1 ~ ln.dist , weights=pp, data = sis,type="BC") summary(gy_logit) sis$trans.overlap1 = (sis$overlap.min.res0.5 * (nrow(sis) - 1) + 0.5) / nrow(sis) gy_logit <- betareg(trans.overlap1 ~ ln.dist , weights=pp,data = sis,type="BC") summary(gy_logit) sis$trans.overlap1 = (sis$overlap.min.res0.1 * (nrow(sis) - 1) + 0.5) / nrow(sis) gy_logit <- betareg(trans.overlap1 ~ ln.dist , weights=pp, data = sis,type="BC") summary(gy_logit) sis$trans.overlap1 = (sis$overlap.min.res0.05 * (nrow(sis) - 1) + 0.5) / nrow(sis) gy_logit <- betareg(trans.overlap1 ~ ln.dist , weights=pp, data = sis,type="BC") summary(gy_logit) ################################################### ### Step 9. Do single genus deletions impact step 7 significance? ################################################### sis$trans.overlap1 = (sis$overlap.min.res1 * (nrow(sis) - 1) + 0.5) / nrow(sis) gens=as.character(unique(sort(sis$genus))) #genus list P1 = c() #make lists for p values for (i in 1:length(gens)) { sis.sub=subset(sis, !(genus==gens[i])) sis.sub$genus=factor(sis.sub$genus) sis.sub$trans.overlap1 = (sis.sub$overlap.min.res1 * (nrow(sis) - 1) + 0.5) / nrow(sis) m1= betareg(trans.overlap1 ~ ln.dist , weights=pp, data=sis.sub, type="BC") sum=summary(m1) my.p1 = round(sum$coefficients$mean[2,4],digits=4) #p value s-o P1=c(P1,my.p1) } mytable = cbind(P1) rownames(mytable) = gens mytable ################################################### ### step 10. Plot relationship between co-occurrence and sister pair age ################################################### sis$trans.overlap1 = (sis$overlap.min.res1 * (nrow(sis) - 1) + 0.5) / nrow(sis) gy_logit <- betareg(trans.overlap1 ~ ln.dist , weights=pp, data = sis,type="BC") summary(gy_logit) tiff("AJB.Fig3.tiff", height = 4, width = 4.5, units = 'in', res=800) par(mar=c(4.3,4.3,1,1), mfrow=c(1,1),cex.lab=1, oma=c(0,0,0,0)) plot(trans.overlap1 ~ ln.dist, data = sis, type = "n", ylab = "Co-occurrence", xlab = "Divergence time (mya)", main = "",xaxt='n', cex.lab=1.5) axis(side=1, at=c(-1.05,0,1.1,2.1,3,3.91), labels=c("0.35","1","3","8","20","50"), cex.axis=0.8) points(trans.overlap1 ~ ln.dist, data = sis, cex = 1, pch = c(19,17,15)[as.numeric(mate.diff)], col = c("gray30","red","pink")[as.numeric(mate.diff)]) points(trans.overlap1 ~ ln.dist, data = sis, cex = 1,pch = c(1,2,0)[as.numeric(mate.diff)]) legend("topright", c("o-o","s-o", "s-s"), title = "", col = c("gray30","red","pink"), pch = c(19,17,15), cex = 1,bty = "y") lines(seq(from = -2, to = 4, by = 0.01), predict(gy_logit, newdata = data.frame(ln.dist=seq(from = -2, to = 4, by = 0.01) )), col = "black", lwd = 2) dev.off() ########################################## ############ Step 11. Does co-occurrence vary by an interaction between sister pair age and mating system? ######## ########################################## sis$trans.overlap1 = (sis$overlap.min.res1 * (nrow(sis) - 1) + 0.5) / nrow(sis) M3 <- betareg(trans.overlap1 ~ mate.diff*ln.dist , weights=pp,data = sis, type="BC") summary(M3) sis$trans.overlap1 = (sis$overlap.min.res0.5 * (nrow(sis) - 1) + 0.5) / nrow(sis) M3 <- betareg(trans.overlap1 ~ mate.diff*ln.dist , weights=pp,data = sis, type="BC") summary(M3) sis$trans.overlap1 = (sis$overlap.min.res0.1 * (nrow(sis) - 1) + 0.5) / nrow(sis) M3 <- betareg(trans.overlap1 ~ mate.diff*ln.dist , weights=pp,data = sis, type="BC") summary(M3) sis$trans.overlap1 = (sis$overlap.min.res0.05 * (nrow(sis) - 1) + 0.5) / nrow(sis) M3 <- betareg(trans.overlap1 ~ mate.diff*ln.dist , weights=pp,data = sis, type="BC") summary(M3) ########################################## ############ Step 12. Histogram of co-occurence for each sister-pair mating system category ########################################## tiff("AJB.Fig2.tiff", height = 4, width = 5, units = 'in', res=800) s.s=subset(sis,(mate.diff=="s_s")) s.o=subset(sis,(mate.diff=="s_o")) o.o=subset(sis,(mate.diff=="o_o")) par(mfrow = c(3,2), oma = c(5,4,3,0) + 0.1, mar = c(1,1,1,3) + 0.1) h = hist(o.o$overlap.min.res0.05, breaks=6, plot=F) h$density = h$counts/sum(h$counts) plot(h, main="", xlab="",ylab='',freq=F, xlim=c(0,1),ylim=c(0,0.8),cex.lab=1.3, axes=F, col="gray30") #axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") axis(side=2,at=c(0,0.2,0.4,0.6,0.8),labels=T, col="darkgray", col.axis="darkgray") abline(v=mean(o.o$overlap.min.res0.05),col="red",lty=2) h = hist(o.o$overlap.min.res1, breaks=15, plot=F) h$density = h$counts/sum(h$counts) plot(h, main="", xlab="",ylab='',freq=F, xlim=c(0,1),ylim=c(0,0.8),cex.lab=1.3, axes=F, col="gray30") #axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") axis(side=2,at=c(0,0.2,0.4,0.6,0.8),labels=T, col="darkgray", col.axis="darkgray") abline(v=mean(o.o$overlap.min.res1),col="red",lty=2) h = hist(s.o$overlap.min.res0.05, breaks=20, plot=F) h$density = h$counts/sum(h$counts) plot(h, main="", xlab="",ylab='',freq=F, xlim=c(0,1),ylim=c(0,0.8),cex.lab=1.3, axes=F, col="red") #axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") axis(side=2,at=c(0,0.2,0.4,0.6,0.8),labels=T, col="darkgray", col.axis="darkgray") abline(v=mean(s.o$overlap.min.res0.05),col="red",lty=2) h = hist(s.o$overlap.min.res1, breaks=15, plot=F) h$density = h$counts/sum(h$counts) plot(h, main="", xlab="",ylab='',freq=F, xlim=c(0,1),ylim=c(0,0.8),cex.lab=1.3, axes=F, col="red") #axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") axis(side=2,at=c(0,0.2,0.4,0.6,0.8),labels=T, col="darkgray", col.axis="darkgray") abline(v=mean(s.o$overlap.min.res1),col="red",lty=2) h = hist(s.s$overlap.min.res0.05, breaks=10, plot=F) h$density = h$counts/sum(h$counts) plot(h, main="", xlab="",ylab='',freq=F, xlim=c(0,1),ylim=c(0,0.8),cex.lab=1.3, axes=F, col="pink") #axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") axis(side=2,at=c(0,0.2,0.4,0.6,0.8),labels=T, col="darkgray", col.axis="darkgray") axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") abline(v=mean(s.s$overlap.min.res0.05),col="red",lty=2) h = hist(s.s$overlap.min.res1, breaks=15, plot=F) h$density = h$counts/sum(h$counts) plot(h, main="", xlab="",ylab='',freq=F, xlim=c(0,1),ylim=c(0,0.8),cex.lab=1.3, axes=F, col="pink") #axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") axis(side=2,at=c(0,0.2,0.4,0.6,0.8),labels=T, col="darkgray", col.axis="darkgray") axis(side=1,at=c(0,0.2,0.4,0.6,0.8,1),labels=T, col="darkgray", col.axis="darkgray") abline(v=mean(s.s$overlap.min.res1),col="red",lty=2) mtext("Co-occurance",side=1,outer=TRUE,line=2.2, cex=1) mtext("Proportion of sister pairs",side=2,outer=TRUE,line=2.5, cex=1) mtext("Finest spatial scale",side=3,outer=TRUE,line=1, cex=1.1, adj=0.1) mtext("(0.05 decimal degree)",side=3,outer=TRUE,line=-0.2, cex=0.7, adj=0.12) mtext("Coarsest spatial scale",side=3,outer=TRUE,line=1, cex=1.1, adj=0.85) mtext("(1 decimal degree)",side=3,outer=TRUE,line=-0.2, cex=0.7, adj=0.8) mtext("o-o",side=4,outer=TRUE,line=-2, cex=1.2, adj=0.82) mtext("s-o",side=4,outer=TRUE,line=-2, cex=1.2, adj=0.47) mtext("s-s",side=4,outer=TRUE,line=-2, cex=1.2, adj=0.1) dev.off()