# Test GCTA maternal effects model with simulated data to using Open Mx ### Testing GCTA approach to maternal effects using highly artificial examples rm(list=ls()) SNPGEN<-function(p,nsnp) { # Function to set up maternal and offspring alleles and genotypes M1<-vector(mode="numeric",length=nsnp) # First maternal allele M2<-vector(mode="numeric",length=nsnp) # Second maternal allele C1<-vector(mode="numeric",length=nsnp) # Maternally derived child allele C2<-vector(mode="numeric",length=nsnp) # Paternally derived child allele # Maternal alleles rand<-runif(nsnp,0,1) M1<-ifelse(rand>p,0,1) rand<-runif(nsnp,0,1) M2<-ifelse(rand>p,0,1) # Paternally-derived child allele rand<-runif(nsnp,0,1) C2<-ifelse(rand>p,0,1) # Select maternally derived child allele rand<-runif(nsnp,0,1) C1<-ifelse(rand>0.5,M1,M2) GM<-M1+M2 GC<-C1+C2 Gtypes<-matrix(nrow=2,ncol=nsnp) Gtypes<-rbind(GM,GC) Gtypes } ################################################################################################ # # # SIMULATE DATA # # # ################################################################################################ # Main program require(OpenMx) set.seed(1932654) # Seeds random number generator (change integer to change sequence) nsnp<-200 # # of SNPS N<-1000 # # of mother-child pairs sde<-2. # se of non-genetic (residual effects) # The snps are divided into those that have direct and indirect effects on offspring, and those that have no effect # Identify samples of SNPs that contribute maternal and offspring effects # Initialize all SNP types GQO<-rep(0,nsnp) GQM<-rep(0,nsnp) noff1<-1 # No. of first SNP contributions to direct effects noff2<-125 # No. of last SNP contributin to direct effects nmom1<-76 # No. of first SNP contribution to indirect ("maternal") effects nmom2<-200 # No. of last SNP contributing to indirect effects. GQO[noff1:noff2]<-1 # Flag loci having offspring effect GQM[nmom1:nmom2]<-1 # Flag loci having maternal effect NO<-sum(GQO) # Observed number of alleles with non-zero offspring effects NM<-sum(GQM) # Observed number of alleles with non-zero maternal effects #GQO+GQM # Check pattern of allelic contributions (1=mother or child, 2=both, 0-null) # Give range of additive effect sizes # Offspring GOlo<-0.45 # Smallest GOhi<-0.55 # Largest # Maternal GMlo<-0.45 # Smallest GMhi<-0.55 # Largest ## Initialize Time to time program sections TIME1 <- proc.time() # Generate additive genetic effects across genome (many of these may subsequently set to zero) GAO<-runif(nsnp,GOlo,GOhi) #Offspring GAM<-runif(nsnp,GMlo,GMhi) #Maternal # Modify effects to exclude loci with no direct or indirect effect GAO<-GAO*GQO # Direct GAM<-GAM*GQM # Indirect # Generate increasing allele frequencies for SNPs p<-runif(nsnp,0.4,0.6) # These are atypically high # Create matrices to hold maternal and offspring genotypes GM<-matrix(nrow=N,ncol=nsnp) # maternal GC<-matrix(nrow=N,ncol=nsnp) # offsping # Create vectors the mother and child raw genetic liabilities PM<-vector(mode="numeric",length=N) # maternal PC<-vector(mode="numeric",length=N) # child Y<-vector(mode="numeric",length=N) # Dyadic phenotype for(pair in (1:N)){ G<-SNPGEN(p,nsnp) GM[pair,1:nsnp]<-G[1,1:nsnp] GC[pair,1:nsnp]<-G[2,1:nsnp] PM[pair]<-sum(G[1,]*GAM) PC[pair]<-sum(G[2,]*GAO) } sdm<-sd(PM) mum<-mean(PM) sdc<-sd(PC) muc<-mean(PC) #PM<-(PM-mum)/sdm #PC<-(PC-muc)/sdc # Generate residuals (normal errors) E<-rnorm(N,0,sde) # Generate dyadic phenotypes Y<-PM+PC+E var(PM) # Variance due to maternal effects (M) var(PC) # Variance due to offsring effects (G) var(E) # Residual variance (E) cov(PM,PC) # Covariance of maternal and offspring effects Q<-2*cov(PM,PC) # G-E Covariance (Q) Q mean(Y) var(Y) # Elapsed time so far TIME2<-proc.time() - TIME1 TIME2 ################################################################################################ # # # Calculate GRMs # # # ################################################################################################ # Time next section TIME1<-proc.time() # Now compute empirical genetic correlations between unrelated mothers (Yang et al.) # Compute Yang et al's "A" matrix for maternal and child genotypes and cross-covariances Alpha<-matrix(nrow=N,ncol=N) # Maternal partition (symmetric) Beta<-matrix(nrow=N,ncol=N) # Offspring partition (symmetric) AMC<-matrix(nrow=N,ncol=N) # Mother-Offfspring partition (square, diagonals +/- 0.5) Delta<-matrix(nrow=N,ncol=N) # Mother-Offfspring partition (sum of AMC and transpose, diagonals +/- 1) for(i in 1:N){ Alpha[i,i]=0 Beta[i,i]=0 for (j in 2:i-1){ Alpha[i,j]=0 Beta[i,j]=0 for(k in 1:nsnp){ Alpha[i,j]<-Alpha[i,j] + (GM[i,k]-2*p[k])*(GM[j,k]-2*p[k])/(2*p[k]*(1-p[k])) Beta[i,j]<-Beta[i,j] + (GC[i,k]-2*p[k])*(GC[j,k]-2*p[k])/(2*p[k]*(1-p[k])) } Alpha[i,j]<-Alpha[i,j]/nsnp Beta[i,j]<-Beta[i,j]/nsnp Alpha[j,i]<-Alpha[i,j] Beta[j,i]<-Beta[i,j] } for(k in 1:nsnp){ Alpha[i,i]<-Alpha[i,i] + ( GM[i,k]^2-(1+2*p[k])*GM[i,k]+2*p[k]^2)/(2*p[k]*(1-p[k])) Beta[i,i]<-Beta[i,i] + ( GC[i,k]^2-(1+2*p[k])*GC[i,k]+2*p[k]^2)/(2*p[k]*(1-p[k])) } Alpha[i,i]<-1+Alpha[i,i]/nsnp Beta[i,i]<-1+Beta[i,i]/nsnp } TIME2<-proc.time() - TIME1 TIME2 TIME1<-proc.time() # Compute cross-related covariances between mothers and children. for(i in 1:N){ for (j in 1:N){ AMC[i,j]<-0 for (k in 1:nsnp){ AMC[i,j]<-AMC[i,j] + (GM[i,k]-2*p[k])*(GC[j,k]-2*p[k])/(2*p[k]*(1-p[k])) } AMC[i,j]<-AMC[i,j]/nsnp } } TIME2<-proc.time() - TIME1 TIME2 TIME1<-proc.time() Delta<-matrix(nrow=N,ncol=N) #Sum of Upper and Lower Triangles of Mother-Offfspring partition (square, diagonals +/- 1) for (i in 1:N){ for (j in i:N){ Delta[i,j]<-AMC[i,j]+AMC[j,i] Delta[j,i]<-Delta[i,j] } } TIME2<-proc.time() - TIME1 TIME2 var(Y) Ybar<-mean(Y) Y<-t(Y) names(Y)<-paste("Y",1:N,sep="") Yobs<-data.frame(Y) selvars<-names(Yobs) ################################################################################################ # # # Fit Unconstrained Model # # # ################################################################################################ # Unconstrained model TIME1<-proc.time() MxTest<-mxModel("GCTA", mxData(observed=Yobs, type="raw"), mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=15,label="m",name="M"), mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=15,label="g",name="G"), mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=3,label="q",name="Q"), mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=4,label="e",name="E"), mxMatrix(type="Iden",nrow=N,ncol=N, name="I"), mxMatrix(type="Full",nrow=1,ncol=N, values=rep(1,N),name="Unit") , mxMatrix(type="Full",nrow=1,ncol=1,free=TRUE,values=Ybar,name="mu"), mxMatrix(type="Full",nrow=N,ncol=N, values=Alpha,name="Alpha"), mxMatrix(type="Full",nrow=N,ncol=N, values=Beta,name="Beta"), mxMatrix(type="Full",nrow=N,ncol=N, values=Delta,name="Delta"), mxAlgebra(expression=G%x%Beta+M%x%Alpha+Q%x%Delta+E%x%I,name="ESigma"), mxAlgebra(expression=mu%x%Unit,name="Ey"), mxFitFunctionML(), mxExpectationNormal(covariance="GCTA.ESigma", means="GCTA.Ey",dimnames=selvars) ) Results<-mxRun(MxTest) TIME2<-proc.time() - TIME1 TIME2 summary(Results)