Created
August 31, 2025 01:10
-
-
Save sashagusev/5fbd954ec7857bd431f1f1c6a73104ef to your computer and use it in GitHub Desktop.
Simulating epistatic phenotypes and estimating their additive heritability
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| library("RColorBrewer") | |
| # --- | |
| # Parameters | |
| # --- | |
| # Samples | |
| N = 2e3 | |
| # Epistatic heritability : | |
| hsq_epistatic = 0.5 | |
| # Allele frequencies to use | |
| # Note MAF==1 will test a uniform distribution | |
| freqs = seq(0.1,1,by=0.1) | |
| # Number of iterations to run | |
| # use 10 for debugging, ~100 for precision | |
| seeds = 10 | |
| # colors | |
| clr = tail(brewer.pal(6,"Reds"),n=3) | |
| # --- | |
| # plot two panels with/without dominance | |
| par(mfrow=c(1,2)) | |
| for(with_dom in c(TRUE,FALSE) ) { | |
| ctr = 1 | |
| for ( M in c(2,5,50) ) { | |
| hsq_joint_add_mean = rep(NA,length(freqs)) | |
| hsq_joint_epi_mean = rep(NA,length(freqs)) | |
| hsq_joint_add_se = rep(NA,length(freqs)) | |
| hsq_joint_epi_se = rep(NA,length(freqs)) | |
| for ( m in 1:length(freqs) ) { | |
| maf = freqs[m] | |
| hsq_joint_add = rep(NA,seeds) | |
| hsq_joint_epi = rep(NA,seeds) | |
| hsq_true = rep(NA,seeds) | |
| for ( s in 1:seeds ) { | |
| if ( maf == 1 ) { | |
| frq = runif(M,0.1,0.9) | |
| X = matrix(nrow=N,ncol=M) | |
| for ( i in 1:ncol(X) ) { | |
| X[,i] = rbinom(N,2,frq[i]) | |
| } | |
| } else { | |
| X = (matrix(rbinom(N*M,2,maf),nrow=N,ncol=M)) | |
| } | |
| # Simulate an epistatic model | |
| gv_epistatic = rep(0,nrow(X)) | |
| for ( i in 1:ncol(X) ) { | |
| bs = rnorm(M) | |
| if ( with_dom ) { | |
| gv_epistatic = gv_epistatic + (X[,i] * X) %*% bs | |
| } else { | |
| gv_epistatic = gv_epistatic + (X[,i] * X) %*% bs - X[,i] * X[,i] * bs[i] | |
| } | |
| } | |
| y = scale(gv_epistatic) * sqrt(hsq_epistatic) + rnorm(N,0,sqrt(1 - hsq_epistatic)) | |
| # True variance explained by the epistatic genetic component | |
| hsq_true[s] = cor(y,gv_epistatic)^2 | |
| # estimate heritability using Haseman-Elston regression | |
| X = scale(X) | |
| k = X %*% t(X) / ncol(X) | |
| yprod = scale( y ) %*% t( scale( y ) ) | |
| yvals = yprod[ lower.tri(yprod,diag=F) ] | |
| kvals = k[ lower.tri(k,diag=F) ] | |
| # Estimated additive heritability | |
| est = lm( yvals ~ kvals ) | |
| hsq_joint_add[s] = est$coef[2] | |
| } | |
| # Store the results | |
| cat( maf , mean(hsq_true) , mean(hsq_joint_add) , sd(hsq_joint_add)/sqrt(seeds) , mean(hsq_joint_epi) , sd(hsq_joint_epi)/sqrt(seeds) , '\n' ) | |
| hsq_joint_add_mean[m] = mean(hsq_joint_add/hsq_true) | |
| hsq_joint_add_se[m] = sd(hsq_joint_add/hsq_true)/sqrt(seeds) | |
| hsq_joint_epi_mean[m] = mean(hsq_joint_epi/hsq_true) | |
| hsq_joint_epi_se[m] = sd(hsq_joint_epi/hsq_true)/sqrt(seeds) | |
| } | |
| # Code for visualization | |
| if(ctr==1) { | |
| if ( with_dom ) plot( 0 , 0 , xaxt="n" , type="n" , xlim=c(0,1) , ylim=c(0,1) , xlab="Allele Frequency" , ylab="Fraction of heritability estimated to be additive" , las=1 , bty="n" , main="Dominance + Epistasis" ) | |
| else plot( 0 , 0 , xaxt="n" , type="n" , xlim=c(0,1) , ylim=c(0,1) , xlab="Allele Frequency" , ylab="Fraction of heritability estimated to be additive" , las=1 , bty="n" , main="Epistasis" ) | |
| axis( side=1 , las=2 , at=freqs , labels = c(freqs[1:9],"Uniform") ) | |
| } | |
| abline(h=1,lty=3) | |
| lines( freqs[1:9] , hsq_joint_add_mean[1:9] , col=clr[ctr] ) | |
| arrows( freqs[1:9] , (hsq_joint_add_mean - hsq_joint_add_se*1.96)[1:9] , freqs[1:9] , (hsq_joint_add_mean + hsq_joint_add_se*1.96)[1:9] , len=0 , col=clr[ctr] ) | |
| points( freqs[1:9] , hsq_joint_add_mean[1:9] , col=clr[ctr] , pch=19 ) | |
| arrows( freqs[10] , (hsq_joint_add_mean - hsq_joint_add_se*1.96)[10] , freqs[10] , (hsq_joint_add_mean + hsq_joint_add_se*1.96)[10] , len=0 , col=clr[ctr] ) | |
| points( freqs[10] , hsq_joint_add_mean[10] , col=clr[ctr] , pch=19 ) | |
| ctr = ctr + 1 | |
| } | |
| legend("bottomright",pch=19,col=clr,legend=c("2","5","50"),bty="n",lwd=1,title="Causal Variants") | |
| } |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment