#=====================================================================================================================================# ### Example MCMCglmm code for VA, VM, IA, and h2 estimation #-------------------------------------------------------------------------------------------------------------------------------------# library(MCMCglmm) data<-read.table("C:/Users/Desktop/R/BC.emer.lgth.txt",na.strings="na",header=T) attach(data) names(data) data$animal<-as.factor(data$animal) data$dam<-as.factor(data$dam) data$sire<-as.factor(data$sire) data$emer.lngth<-as.numeric(data$emer.lgth) head(data) ped<-read.table("C:/Users/Desktop/R/BC.emer.pedigree.txt",na.strings="na",header=T) attach(ped) names(ped) head(ped) prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped,data=data,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h2<-model$VCV[,"animal"]/(model$VCV[,"animal"]+model$VCV[,"dam"]+model$VCV[,"units"]) mean.emer<-mean(emer.lgth,na.rm=TRUE) posterior.IA<-(model$VCV[,"animal"])/(mean.emer^2) #EMERGENCE LENGTH plot(model$Sol) plot(model$VCV) plot(posterior.h2) plot(posterior.IA) autocorr(model$VCV) effectiveSize(model$Sol) effectiveSize(model$VCV) effectiveSize(posterior.h2) effectiveSize(posterior.IA) heidel.diag(model$VCV) posterior.mode(model$Sol) HPDinterval(model$Sol,0.95) posterior.mode(posterior.h2) HPDinterval(posterior.h2,0.95) posterior.mode(posterior.IA) HPDinterval(posterior.IA,0.95) posterior.mode(model$VCV) HPDinterval(model$VCV) summary(model) #=====================================================================================================================================# ### Example MCMCglmm code for VM estimation for survival #-------------------------------------------------------------------------------------------------------------------------------------# library(MCMCglmm) data<-read.table("C:/Users/Desktop/R/BC.survival.txt",na.strings="na",header=T) attach(data) names(data) data$animal<-as.factor(data$animal) data$dam<-as.factor(data$dam) data$sire<-as.factor(data$sire) data$survived<-as.factor(data$survived) ped<-read.table("C:/Users/Jacquelyn/Desktop/R/Chapter 3/BC.surv.pedigree.txt",na.strings="na",header=T) attach(ped) names(ped) head(ped) prior = list(R=list(V=1,nu=0,fix=1),G=list(G1=list(V=0.5,nu=0.0164),G2=list(V=0.5,nu=0.0164))) model<-MCMCglmm(survived~1,random=~dam+dam:sire,family="categorical",pedigree=ped,rcov=~units,data=data,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=T) plot(model$Sol) plot(model$VCV) autocorr(model$VCV) summary(model) posterior.mode(model$VCV) HPDinterval(model$VCV) #=====================================================================================================================================# ### Example code for estimating the posterior distribution of r #-------------------------------------------------------------------------------------------------------------------------------------# library(MCMCglmm) #BC data1<-read.table("C:/Users/Desktop/R/BC.emer.txt",na.strings="na",header=T) attach(data1) names(data1) data1$animal<-as.factor(data1$animal) data1$dam<-as.factor(data1$dam) data1$sire<-as.factor(data1$sire) data1$emer.lgth<-as.numeric(data1$emer.lgth) head(data1) ped1<-read.table("C:/Users/Desktop/R/BC.emer.pedigree.txt",na.strings="na",header=T) attach(ped1) names(ped1) head(ped1) #-------------------------------------------------------------- #BC EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model1<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped1,data=data1,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h21<-model1$VCV[,"animal"]/(model1$VCV[,"animal"]+model1$VCV[,"dam"]+model1$VCV[,"units"]) mean.emer1<-mean(emer.lgth,na.rm=TRUE) posterior.IA1<-(model1$VCV[,"animal"])/(mean.emer1^2) #-------------------------------------------------------------- #BF data2<-read.table("C:/Users/Desktop/R/BF.emer.txt",na.strings="na",header=T) attach(data2) names(data2) data2$animal<-as.factor(data2$animal) data2$dam<-as.factor(data2$dam) data2$sire<-as.factor(data2$sire) data2$emer.lgth<-as.numeric(data2$emer.lgth) head(data2) ped2<-read.table("C:/Users/Desktop/R/BF.emer.pedigree.txt",na.strings="na",header=T) attach(ped2) names(ped2) head(ped2) #-------------------------------------------------------------- #BF EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model2<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped2,data=data2,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h22<-model2$VCV[,"animal"]/(model2$VCV[,"animal"]+model2$VCV[,"dam"]+model2$VCV[,"units"]) mean.emer2<-mean(emer.lgth,na.rm=TRUE) posterior.IA2<-(model2$VCV[,"animal"])/(mean.emer2^2) #-------------------------------------------------------------- #CC data3<-read.table("C:/Users/Desktop/R/CC.emer.txt",na.strings="na",header=T) attach(data3) names(data3) data3$animal<-as.factor(data3$animal) data3$dam<-as.factor(data3$dam) data3$sire<-as.factor(data3$sire) data3$emer.lgth<-as.numeric(data3$emer.lgth) head(data3) ped3<-read.table("C:/Users/Desktop/R/CC.emer.pedigree.txt",na.strings="na",header=T) attach(ped3) names(ped3) head(ped3) #-------------------------------------------------------------- #CC EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model3<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped3,data=data3,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h23<-model3$VCV[,"animal"]/(model3$VCV[,"animal"]+model3$VCV[,"dam"]+model3$VCV[,"units"]) mean.emer3<-mean(emer.lgth,na.rm=TRUE) posterior.IA3<-(model3$VCV[,"animal"])/(mean.emer3^2) #------------------------------------------------------------- #DY data4<-read.table("C:/Users/Desktop/R/DY.emer.txt",na.strings="na",header=T) attach(data4) names(data4) data4$animal<-as.factor(data4$animal) data4$dam<-as.factor(data4$dam) data4$sire<-as.factor(data4$sire) data4$emer.lgth<-as.numeric(data4$emer.lgth) head(data4) ped4<-read.table("C:/Users/Desktop/R/DY.emer.pedigree.txt",na.strings="na",header=T) attach(ped4) names(ped4) head(ped4) #-------------------------------------------------------------- #DY EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model4<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped4,data=data4,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h24<-model4$VCV[,"animal"]/(model4$VCV[,"animal"]+model4$VCV[,"dam"]+model4$VCV[,"units"]) mean.emer4<-mean(emer.lgth,na.rm=TRUE) posterior.IA4<-(model4$VCV[,"animal"])/(mean.emer4^2) #-------------------------------------------------------------- #FW data5<-read.table("C:/Users/Desktop/R/FW.emer.txt",na.strings="na",header=T) attach(data5) names(data5) data5$animal<-as.factor(data5$animal) data5$dam<-as.factor(data5$dam) data5$sire<-as.factor(data5$sire) data5$emer.lgth<-as.numeric(data5$emer.lgth) head(data5) ped5<-read.table("C:/Users/Desktop/R/FW.emer.pedigree.txt",na.strings="na",header=T) attach(ped5) names(ped5) head(ped5) #-------------------------------------------------------------- #FW EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model5<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped5,data=data5,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h25<-model5$VCV[,"animal"]/(model5$VCV[,"animal"]+model5$VCV[,"dam"]+model5$VCV[,"units"]) mean.emer5<-mean(emer.lgth,na.rm=TRUE) posterior.IA5<-(model5$VCV[,"animal"])/(mean.emer5^2) #-------------------------------------------------------------- #STBC data6<-read.table("C:/Users/Desktop/R/STBC.emer.txt",na.strings="na",header=T) attach(data6) names(data6) data6$animal<-as.factor(data6$animal) data6$dam<-as.factor(data6$dam) data6$sire<-as.factor(data6$sire) data6$emer.lgth<-as.numeric(data6$emer.lgth) head(data6) ped6<-read.table("C:/Users/Desktop/R/STBC.emer.pedigree.txt",na.strings="na",header=T) attach(ped6) names(ped6) head(ped6) #-------------------------------------------------------------- #STBC EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model6<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped6,data=data6,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h26<-model6$VCV[,"animal"]/(model6$VCV[,"animal"]+model6$VCV[,"dam"]+model6$VCV[,"units"]) mean.emer6<-mean(emer.lgth,na.rm=TRUE) posterior.IA6<-(model6$VCV[,"animal"])/(mean.emer6^2) #-------------------------------------------------------------- #UO data7<-read.table("C:/Users/Desktop/R/UO.emer.txt",na.strings="na",header=T) attach(data7) names(data7) data7$animal<-as.factor(data7$animal) data7$dam<-as.factor(data7$dam) data7$sire<-as.factor(data7$sire) data7$emer.lgth<-as.numeric(data7$emer.lgth) head(data7) ped7<-read.table("C:/Users/Desktop/R/UO.emer.pedigree.txt",na.strings="na",header=T) attach(ped7) names(ped7) head(ped7) #-------------------------------------------------------------- #UO EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model7<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped7,data=data7,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h27<-model7$VCV[,"animal"]/(model7$VCV[,"animal"]+model7$VCV[,"dam"]+model7$VCV[,"units"]) mean.emer7<-mean(emer.lgth,na.rm=TRUE) posterior.IA7<-(model7$VCV[,"animal"])/(mean.emer7^2) #-------------------------------------------------------------- #WC data8<-read.table("C:/Users/Desktop/R/WC.emer.txt",na.strings="na",header=T) attach(data8) names(data8) data8$animal<-as.factor(data8$animal) data8$dam<-as.factor(data8$dam) data8$sire<-as.factor(data8$sire) data8$emer.lgth<-as.numeric(data8$emer.lgth) head(data8) ped8<-read.table("C:/Users/Desktop/R/WC.emer.pedigree.txt",na.strings="na",header=T) attach(ped8) names(ped8) head(ped8) #-------------------------------------------------------------- #WC EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model8<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped8,data=data8,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h28<-model8$VCV[,"animal"]/(model8$VCV[,"animal"]+model8$VCV[,"dam"]+model8$VCV[,"units"]) mean.emer8<-mean(emer.lgth,na.rm=TRUE) posterior.IA8<-(model8$VCV[,"animal"])/(mean.emer8^2) #-------------------------------------------------------------- #WN data9<-read.table("C:/Users/Desktop/R/WN.emer.txt",na.strings="na",header=T) attach(data9) names(data9) data9$animal<-as.factor(data9$animal) data9$dam<-as.factor(data9$dam) data9$sire<-as.factor(data9$sire) data9$emer.lgth<-as.numeric(data9$emer.lgth) head(data9) ped9<-read.table("C:/Users/Desktop/R/WN.emer.pedigree.txt",na.strings="na",header=T) attach(ped9) names(ped9) head(ped9) #-------------------------------------------------------------- #WN EMERGENCE LENGTH prior<-list(G=list(G1=list(V=1,nu=1),G2=list(V=1,nu=1)),R=list(V=1,nu=1)) model9<-MCMCglmm(emer.lgth~1,random=~animal+dam,pedigree=ped9,data=data9,nitt=1000000,thin=50,burnin=300000,prior=prior,verbose=FALSE) posterior.h29<-model9$VCV[,"animal"]/(model9$VCV[,"animal"]+model9$VCV[,"dam"]+model9$VCV[,"units"]) mean.emer9<-mean(emer.lgth,na.rm=TRUE) posterior.IA9<-(model9$VCV[,"animal"])/(mean.emer9^2) #--------------------------------------------------------------------- EMERGENCE LENGTH CORRELATIONS Post.VA <-cbind(model1$VCV[,"animal"],model2$VCV[,"animal"],model3$VCV[,"animal"],model4$VCV[,"animal"],model5$VCV[,"animal"],model6$VCV[,"animal"],model7$VCV[,"animal"],model8$VCV[,"animal"],model9$VCV[,"animal"]) Post.h2 <-cbind(posterior.h21,posterior.h22,posterior.h23,posterior.h24,posterior.h25,posterior.h26,posterior.h27,posterior.h28,posterior.h29) Post.IA <-cbind(posterior.IA1,posterior.IA2,posterior.IA3,posterior.IA4,posterior.IA5,posterior.IA6,posterior.IA7,posterior.IA8,posterior.AI9) Post.VM <-cbind(model1$VCV[,"dam"],model2$VCV[,"dam"],model3$VCV[,"dam"],model4$VCV[,"dam"],model5$VCV[,"dam"],model6$VCV[,"dam"],model7$VCV[,"dam"],model8$VCV[,"dam"],model9$VCV[,"dam"]) N.table<-read.table("C:/Users/Desktop/R/N.txt",na.strings="na",header=T) attach(N.table) names(N.table) Nb.table<-read.table("C:/Users/Desktop/R/Nb.txt",na.strings="na",header=T) attach(Nb.table) names(Nb.table) A <- as.matrix(Post.VA) B <- as.matrix(N.table) C <- as.matrix(Nb.table) D <- as.matrix(Post.h2) E <- as.matrix(Post.IA) F <- as.matrix(Post.VM) Corr.resultsVA.N<-sapply(seq.int(dim(A)[1]), function(i) cor(A[i,], B[i,])) Corr.resultsVA.Nb<-sapply(seq.int(dim(A)[1]), function(i) cor(A[i,], C[i,])) Corr.resultsh2.N<-sapply(seq.int(dim(D)[1]), function(i) cor(D[i,], B[i,])) Corr.resultsh2.Nb<-sapply(seq.int(dim(D)[1]), function(i) cor(D[i,], C[i,])) Corr.resultsIA.N<-sapply(seq.int(dim(E)[1]), function(i) cor(E[i,], B[i,])) Corr.resultsIA.Nb<-sapply(seq.int(dim(E)[1]), function(i) cor(E[i,], C[i,])) Corr.resultsVM.N<-sapply(seq.int(dim(F)[1]), function(i) cor(F[i,], B[i,])) Corr.resultsVM.Nb<-sapply(seq.int(dim(F)[1]), function(i) cor(F[i,], C[i,])) ##VA Corr.resultsVA.N hist(Corr.resultsVA.N) mcmc(data= Corr.resultsVA.N, start = 1, end = 14000, thin = 1) Corr.resultsVA.N.mcmc<-as.mcmc(Corr.resultsVA.N) HPDinterval(Corr.resultsVA.N.mcmc, 0.95) posterior.mode(Corr.resultsVA.N.mcmc) Corr.resultsVA.Nb hist(Corr.resultsVA.Nb) mcmc(data= Corr.resultsVA.Nb, start = 1, end = 14000, thin = 1) Corr.resultsVA.Nb.mcmc<-as.mcmc(Corr.resultsVA.Nb) HPDinterval(Corr.resultsVA.Nb.mcmc, 0.95) posterior.mode(Corr.resultsVA.Nb.mcmc) ##h2 Corr.resultsh2.N hist(Corr.resultsh2.N) mcmc(data= Corr.resultsh2.N, start = 1, end = 14000, thin = 1) Corr.resultsh2.N.mcmc<-as.mcmc(Corr.resultsh2.N) HPDinterval(Corr.resultsh2.N.mcmc, 0.95) posterior.mode(Corr.resultsh2.N.mcmc) Corr.resultsh2.Nb hist(Corr.resultsh2.Nb) mcmc(data= Corr.resultsh2.Nb, start = 1, end = 14000, thin = 1) Corr.resultsh2.Nb.mcmc<-as.mcmc(Corr.resultsh2.Nb) HPDinterval(Corr.resultsh2.Nb.mcmc, 0.95) posterior.mode(Corr.resultsh2.Nb.mcmc) ##IA Corr.resultsIA.N hist(Corr.resultsIA.N) mcmc(data= Corr.resultsIA.N, start = 1, end = 14000, thin = 1) Corr.resultsIA.N.mcmc<-as.mcmc(Corr.resultsIA.N) HPDinterval(Corr.resultsIA.N.mcmc, 0.95) posterior.mode(Corr.resultsIA.N.mcmc) Corr.resultsIA.Nb hist(Corr.resultsIA.Nb) mcmc(data= Corr.resultsIA.Nb, start = 1, end = 14000, thin = 1) Corr.resultsIA.Nb.mcmc<-as.mcmc(Corr.resultsIA.Nb) HPDinterval(Corr.resultsIA.Nb.mcmc, 0.95) posterior.mode(Corr.resultsIA.Nb.mcmc) ##VM Corr.resultsVM.N hist(Corr.resultsVM.N) mcmc(data= Corr.resultsVM.N, start = 1, end = 14000, thin = 1) Corr.resultsVM.N.mcmc<-as.mcmc(Corr.resultsVM.N) HPDinterval(Corr.resultsVM.N.mcmc, 0.95) posterior.mode(Corr.resultsVM.N.mcmc) Corr.resultsVM.Nb hist(Corr.resultsVM.Nb) mcmc(data= Corr.resultsVM.Nb, start = 1, end = 14000, thin = 1) Corr.resultsVM.Nb.mcmc<-as.mcmc(Corr.resultsVM.Nb) HPDinterval(Corr.resultsVM.Nb.mcmc, 0.95) posterior.mode(Corr.resultsVM.Nb.mcmc) #=====================================================================================================================================# ### Example MCMCglmm code for QST estimation #-------------------------------------------------------------------------------------------------------------------------------------# library(MCMCglmm) dat<-read.table("C:/Users/Desktop/R/BCxBF.emer.txt",na.strings="na",header=T) attach(dat) names(dat) dat$population<-as.factor(dat$population) dat$dam<-as.factor(dat$dam) dat$sire<-as.factor(dat$sire) dat$emer.lgth<-as.numeric(dat$emer.lgth) names(dat)<-c("Stream","Pop","DAM","SIRE","emer.lgth") dat$DAM<-as.factor(paste(rep("F", dim(dat)[1]), dat$DAM, sep=".")) dat$SIRE<-as.factor(paste(rep("M", dim(dat)[1]), dat$SIRE, sep=".")) dat$animal<-as.factor(paste(dat$Pop, 1:dim(dat)[1], sep=".")) ped<-data.frame(ID=c(levels(dat$DAM), levels(dat$SIRE))) ped$DAM<-rep(NA, dim(ped)[1]) ped$SIRE<-rep(NA, dim(ped)[1]) ped<-rbind(ped, data.frame(ID=dat$animal, DAM=dat$DAM, SIRE=dat$SIRE)) p.var<-var(dat$emer.lgth,na.rm=T) prior.full<-list(G=list(G1=list(V=matrix(p.var/4),n=1), G2=list(V=matrix(p.var/4),n=1), G3=list(V=matrix(p.var/4),n=1)),R=list(V=matrix(p.var/4),n=1)) model.full<-MCMCglmm(emer.lgth~1, random=~animal+DAM+Pop, pedigree=ped, prior=prior.full, data=dat, nitt=1000000, thin=500, burnin=300000,verbose=FALSE) plot(model.full) autocorr.plot(model.full$VCV) autocorr(model.full$VCV) plot(model.full$VCV[,"Pop"]) autocorr.plot(model.full$VCV[,"Pop"]) posterior.mode(model.full$VCV[,"Pop"]) HPDinterval(model.full$VCV[,"Pop"]) plot(model.full$VCV[,"DAM"]) autocorr.plot(model.full$VCV[,"DAM"]) posterior.mode(model.full$VCV[,"DAM"]) HPDinterval(model.full$VCV[,"DAM"]) plot(model.full$VCV[,"animal"]) autocorr.plot(model.full$VCV[,"animal"]) posterior.mode(model.full$VCV[,"animal"]) HPDinterval(model.full$VCV[,"animal"]) qst<-model.full$VCV[,"Pop"]/(model.full$VCV[,"Pop"] + 2*model.full$VCV[,"animal"]) plot(qst) autocorr.plot(qst) posterior.mode(qst) HPDinterval(qst) #=====================================================================================================================================# ### Example code for Mantel and partial Mantel tests #-------------------------------------------------------------------------------------------------------------------------------------# library(vegan) Mean.Nb <- as.matrix(read.table("clipboard", head=T, row.names=1)) Mean.N <- as.matrix(read.table("clipboard", head=T, row.names=1)) Fst.micro <- as.matrix(read.table("clipboard", head=T, row.names=1)) Fst.SNPs <- as.matrix(read.table("clipboard", head=T, row.names=1)) Qst.emer<- as.matrix(read.table("clipboard", head=T, row.names=1)) mantel(as.dist(Mean.Nb), as.dist(Fst.micro), permutations = 10000) mantel(as.dist(Mean.N), as.dist(Fst.micro), permutations = 10000) mantel(as.dist(Mean.Nb), as.dist(Fst.SNPs), permutations = 10000) mantel(as.dist(Mean.N), as.dist(Fst.SNPs), permutations = 10000) mantel.partial(as.dist(Qst.emer), as.dist(Mean.N), as.dist(Fst.micro), permutations = 9999) mantel.partial(as.dist(Qst.emer), as.dist(Mean.N), as.dist(Fst.SNPs), permutations = 9999) mantel.partial(as.dist(Qst.emer), as.dist(Mean.Nb), as.dist(Fst.micro), permutations = 9999) mantel.partial(as.dist(Qst.emer), as.dist(Mean.Nb), as.dist(Fst.SNPs), permutations = 9999)