#Script for performing simulations pin<-function (object, transform) { pframe <- as.list(object$gammas) names(pframe) <- paste("V", seq(1, length(pframe)), sep = "") tvalue <- eval(deriv(transform[[length(transform)]], names(pframe)), pframe) X <- as.vector(attr(tvalue, "gradient")) X[object$gammas.type == 1] <- 0 tname <- if (length(transform) == 3) transform[[2]] else "" n <- length(pframe) i <- rep(1:n, 1:n) j <- sequence(1:n) k <- 1 + (i > j) Vmat <- object$ai se <- sqrt(sum(Vmat * X[i] * X[j] * k)) data.frame(row.names = tname, Estimate = tvalue, SE = se) } #いいいいいいいいいいいいいい ####Start of the script###### #いいいいいいいいいいいいいい Exp_age<-read.table("Pheno_Data_AGE_ESM.txt",header=T) library(asreml) ped<-read.table("Pedigree_AGE_ESM.txt", header=T) ainv= asreml.Ainverse(ped[,c(1,3,2)])$ginv #Open pedantics library(pedantics) Pval.GXA<-data.frame("iteration"=NA, "Converge1"=NA,"Converge2"=NA, "pvalGXA"=NA, "VA"=NA, "VR1"=NA,"VR2"=NA,"VR3"=NA,"VR4"=NA, "Converge3"=NA,"rG"=NA,"SE_rG"=NA,"h2J"=NA,"SE_h2J"=NA,"h2A"=NA,"SE_h2A"=NA, "pvalrG0"=NA,"pvalVAR"=NA) for (i in 1:1000){ Pval.GXA[i,"iteration"]<-i simPhen<-phensim(ped,traits=1, randomA=8,randomE=41)$phenotypes simPhen<-droplevels(simPhen[simPhen$id %in%Exp_age$Ring,]) names(simPhen)[1]<-"Ring" simPhen<-merge(simPhen,Exp_age) #Define age as blocks: 0,1,2,>=3 years old simPhen[simPhen$TESTAGE< 12,"block"]<-1 simPhen[simPhen$TESTAGE>=12 &simPhen$TESTAGE<24 ,"block"]<-2 simPhen[simPhen$TESTAGE>=24 &simPhen$TESTAGE<36 ,"block"]<-3 simPhen[simPhen$TESTAGE>=36,"block"]<-4 #Mean-center age simPhen$AGE<-simPhen$block-mean(simPhen$block,na.rm=T) #Average trait responds to age simPhen$trait_1b<-simPhen$trait_1+(1.87E-01 *simPhen$TESTAGE -2.88E-03*simPhen$TESTAGE^2) #Transform to get a similar distribution of raw data simPhen$trait_1c<-(simPhen$trait_1b^1.1) #simPhen$trait_1c<-simPhen$trait_1c/var(simPhen$trait_1c,na.rm=T) *49 #Separate juveniles and adults simPhen$LOCOMYOUNG<-ifelse(simPhen$TESTAGE <12& simPhen$June.Days <283,simPhen$trait_1b,NA ) simPhen$LOCOMADULT<-ifelse(simPhen$TESTAGE %in% c(12:98),simPhen$trait_1b,NA) #Each individual needs to have one line for each block, even if not measured in each (NAs) blocks<-merge(simPhen$Ring,unique(simPhen$block)) names(blocks)<-c("Ring", "block") simPhen2<-merge(simPhen,blocks,all=T) simPhen2<-simPhen2[order(simPhen2$block),]#Data needs to be sorted by block ##Run the models simPhen2$block<-as.factor(simPhen2$block) simPhen2$Ring<-as.factor(as.character(simPhen2$Ring)) #Test if mean age effects are correct #model.LOCOM <- asreml(trait_1c ~ 1 + TESTAGE + I(TESTAGE^2), # random = ~ped(Ring), # data = simPhen,na.method.X="include", # ginverse = list(Ring=ainv), # maxiter=500) #summary(model.LOCOM)$varcomp #model.LOCOM$coefficients$fixed #wald.asreml(model.LOCOM) ################################################## #Simple model with homegenous residuals and years ################################################## #model.LOCOM00 <- asreml(trait_1c ~ 1 + as.factor(AGE) , # random = ~ped(Ring), # data = simPhen2,na.method.X="include", # ginverse = list(Ring=ainv), # maxiter=100) #summary(model.LOCOM00)$varcomp ################################################## #Simple model with heterogenous residuals and years ################################################## model.LOCOM0a <- asreml(trait_1b ~ 1 + as.factor(AGE) , random = ~ped(Ring), rcov= ~ at(block):units, data = simPhen2,na.method.X="include", ginverse = list(Ring=ainv), maxiter=100) #summary(model.LOCOM0a)$varcomp #1-pchisq(2*abs(model.LOCOM0a$loglik-model.LOCOM00$loglik),3)#Test for heterogeneous vs.homogeneous residuals Pval.GXA[i,"Converge1"]<-model.LOCOM0a$converge Pval.GXA[i,"VA"]<-summary(model.LOCOM0a)$varcomp[1,1] Pval.GXA[i,"VR1"]<-summary(model.LOCOM0a)$varcomp[2,1] Pval.GXA[i,"VR2"]<-summary(model.LOCOM0a)$varcomp[3,1] Pval.GXA[i,"VR3"]<-summary(model.LOCOM0a)$varcomp[4,1] Pval.GXA[i,"VR4"]<-summary(model.LOCOM0a)$varcomp[5,1] ################################################################### #GxA: Model where VA varies as a linear function of age (x=1) ################################################################### model.LOCOM1 <- asreml(trait_1b ~ 1 + as.factor(AGE), random = ~us(pol(AGE,1),init=c(4,0.1,0.1)):ped(Ring), rcov= ~ at(block,init=c(40,40,40,40)):units, data = simPhen2,na.method.X="include", ginverse = list(Ring=ainv), maxiter=100) #summary(model.LOCOM1)$varcomp #model.LOCOM1<-update(model.LOCOM1) Pval.GXA[i,"Converge2"]<-model.LOCOM1$converge Pval.GXA[i,"pvalGXA"]<-ifelse(model.LOCOM1$converge==TRUE,1-pchisq(2*abs(model.LOCOM0a$loglik-model.LOCOM1$loglik),2),NA)#Test for GxA ################################################################### #Character State approach ################################################################### simPhen$Ring<-as.factor(as.character(simPhen$Ring)) model.biv<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + TESTAGE + I(TESTAGE^2), random = ~us(trait,init=c(8,8,10)):ped(Ring), rcov= ~ units:diag(trait,init=c(41,41)), data = simPhen,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=100) Pval.GXA[i,"Converge3"]<-model.biv$converge #summary(model.biv)$varcomp #summary(model.biv.diag) if(model.biv$converge==TRUE){ Pval.GXA[i,"rG"]<-pin(model.biv, rA~V2/sqrt(V1*V3))[1]#Genetic correlation across ages Pval.GXA[i,"h2J"]<-pin(model.biv, h2j~V1/(V1+V5))[1]#Heritability in young Pval.GXA[i,"h2A"]<-pin(model.biv, h2a~V3/(V3+V6))[1]#Heritability in adults Pval.GXA[i,"SE_rG"]<-pin(model.biv, rA~V2/sqrt(V1*V3))[2]#SE of genetic correlation across ages Pval.GXA[i,"SE_h2J"]<-pin(model.biv, h2j~V1/(V1+V5))[2]#SE of heritability in young Pval.GXA[i,"SE_h2A"]<-pin(model.biv, h2a~V3/(V3+V6))[2]#SE of heritability in adults #Test if correlation is different from zero model.biv.diag<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + TESTAGE + I(TESTAGE^2), random = ~diag(trait,init=c(8,8)):ped(Ring), rcov= ~ units:diag(trait,init=c(41,41)), data = simPhen,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=100) Pval.GXA[i,"pvalrG0"]<-1-pchisq(2*abs(model.biv$loglik - model.biv.diag$loglik),1) } #Test if additive genetic variances are the same between young and adults sv<- asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + TESTAGE + I(TESTAGE^2), random = ~diag(trait,init=c(8,8)):ped(Ring), rcov= ~ units:diag(trait,init=c(41,41)), data = simPhen,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=100,start.values = T) gam<-sv$gammas.table gam$fac<-factor(c(1,1,3,4,5))#same VA M <- asreml.constraints(~ fac, gammas=gam) model.biv.VA<-asreml(cbind(LOCOMYOUNG,LOCOMADULT) ~ trait + TESTAGE + I(TESTAGE^2), random = ~diag(trait):ped(Ring), rcov= ~ units:diag(trait), data = simPhen,na.method.X="include",na.method.Y="include", ginverse = list(Ring=ainv), maxiter=100,constraints = M,G.param=gam,R.param=gam) #summary(model.biv.VA) Pval.GXA[i,"pvalVAR"]<-1-pchisq(2*abs(model.biv.diag$loglik - model.biv.VA$loglik),1)# Test if VA(or VR) are statistically different across age groups if( i %in% c(300,400,500,600,700,800,900,1000)){write.table(Pval.GXA,"Results_Sim_No_GXA.txt",row.names = F)}#Save the results regularly } #######End of the sumulation########## #いいいいいいいいいいいいいい #Analyze the results #いいいいいいいいいいいいいい summary(Pval.GXA) #Calculate rate of detection of GxA (GxA does not occur in simulated data) Pval.GXA.red1<-Pval.GXA[Pval.GXA$Converge2==TRUE, ] nrow(Pval.GXA.red1[Pval.GXA.red1$pvalGXA<0.05,])/nrow(Pval.GXA.red1)#5% false positive #Calculate rate of detection of rG across ages (this correlation should be 1) Pval.GXA.red2<-Pval.GXA[Pval.GXA$Converge3==TRUE, ] nrow(Pval.GXA.red2[Pval.GXA.red2$pvalrG0<0.05,])/nrow(Pval.GXA.red2)#62% power to detect a positive genetic correlation across ages #probabilities combined Pval.GXA.red<-Pval.GXA[Pval.GXA$Converge2==TRUE & Pval.GXA$Converge3==TRUE, ] nrow(Pval.GXA.red[Pval.GXA.red$pvalrG0>0.05 & Pval.GXA.red$pvalGXA<0.05,])/nrow(Pval.GXA.red)#4% power to detect no genetic correlation across ages but GxA #Calculate rate of detection of different VAs across ages nrow(Pval.GXA.red2[Pval.GXA.red2$pvalVAR<0.05,])/nrow(Pval.GXA.red2)#3% mean(Pval.GXA.red2$rG)#0.78 mean(Pval.GXA.red2$h2A)#0.21 quantile(Pval.GXA.red2$rG,probs=c(0.025,0.975))#95%CI quantile(Pval.GXA.red2$h2A,probs=c(0.025,0.975))#95%CI sum(Pval.GXA.red2$rG-1.96*(Pval.GXA.red2$SE_rG)>0)/length(Pval.GXA.red2$rG)#rG statistically >0 52% times sum(Pval.GXA.red2$h2A-1.96*(Pval.GXA.red2$SE_h2A)>0)/length(Pval.GXA.red2$h2A)#h2Adults statistically >0 32% times sum(Pval.GXA.red2$h2J-1.96*(Pval.GXA.red2$SE_h2J)>0)/length(Pval.GXA.red2$h2J)#h2Juv statistically >0 99% times