PMR Null Simulation with Sparsity

Overview

This code evaluates PMR under the null when horizontal pleiotropy is assumed to be sparse. In the PMR paper, the authors claim that the false false positive rate is controlled when horizontal pleiotropic effects are sparse (Supplementary Figure 6). In their main simulation, all variants have exactly the same horizontal pleotropic effect (\(\gamma = 0.001\)). In the supplement, they vary the effect size (0.0001, 0.0005, 0.001, 0.002) and assume that only fraction of variants are non-zero (10%, 30%, 50%, 100%). They provide some example code for their simulations and running PMR here. Their simulations hold effect size fixed and then reduce the number of variants to induce sparsity. This means that the PVE of pleiotropy is reduced in the sparse scenarios.

Our cTWAS simulations suggest that the false positive rate of PMR is not controlled. Our simulations assume SNP effects that are much sparser than the 10% in the PMR paper. I also think that the reduced PVE for sparse scenarios in their simulations masks potential issues with the false positive rate. To investigate, I adapted the code provided with the PMR paper to hold PVE fixed when the sparsity is changed, and allowed different numbers of non-zero pleiotropy effects. I also investigated the sparsity of the simulated eQTL effects, though it seemed to not make much difference in false positive rate.

Simulations

library(PMR)

sparse_values <- c(50, 5, 1)
PVE_values <- c(0.0001, 0.001, 0.002, 0.01)

####################

results <- data.frame(sparsity=as.integer(rep(NA, length(sparse_values)*length(PVE_values))),
                      pve=as.integer(rep(NA, length(sparse_values)*length(PVE_values))),
                      gamma=as.integer(rep(NA, length(sparse_values)*length(PVE_values))),
                      gamma_hat=as.integer(rep(NA, length(sparse_values)*length(PVE_values))),
                      false_positive_rate=as.integer(rep(NA, length(sparse_values)*length(PVE_values))))

#load the scaled genenotype matrix in eQTL data (e.g. cis-SNPs of BACE1 gene from GEUVADIS data)
zx<-read.table("/home/wcrouse/PMR/PMR/example/zx.txt")
zx<-as.matrix(zx)
n1 = dim(zx)[1]
q = dim(zx)[2]

#load the scaled genenotype matrix in GWAS data (e.g. the same cis-SNPs from GERA data)
zy<-read.table("/home/wcrouse/PMR/PMR/example/zy.txt")
zy<-as.matrix(zy)
n2<-dim(zy)[1]

#number of iterations for each condition
nsims <- 1000

iter <- 0

for (n_effects_zy in sparse_values){
  for (PVE_zy in PVE_values){
    
    iter <- iter + 1
    
    # print(iter)
    
    #set number of genotype effects in gene expression
    n_effects_zx <- q
    
    #set PVE by genotypes for gene expression
    PVE_zx <- 0.1
    
    alpha_hat <- rep(NA, nsims)
    gamma_hat <- rep(NA, nsims)
    pvalue_alpha <- rep(NA, nsims)
    pvalue_gamma <- rep(NA, nsims)
    Gamma_all <- rep(NA, nsims)
    
    for (i in 1:nsims){
      # if (i %% 10 == 0){
      #   print(i)
      # }
      
      #get the simulated gene expression data
      beta <- matrix(0, q, 1)
      beta[sample(1:length(beta), n_effects_zx), 1] <- rnorm(n_effects_zx)
      
      x_mean <- as.vector(zx%*%beta)
      scaling_factor <- sqrt(PVE_zx)/sd(x_mean)
      x_mean <- x_mean*scaling_factor
      beta <- beta*scaling_factor #compute scaled beta
      
      epsilon_x <- rnorm(n1)
      epsilon_x <- epsilon_x/sd(epsilon_x)*sqrt(1-PVE_zx)
      
      x <- as.matrix(x_mean + epsilon_x)
      
      #get the simulated phenotype
      Gamma <- rep(0, q)
      Gamma[sample(1:length(Gamma), n_effects_zy)] <- 1
      
      y_mean <- as.vector(zy%*%Gamma)
      scaling_factor <- sqrt(PVE_zy)/sd(y_mean)
      y_mean <- y_mean*scaling_factor
      Gamma <- Gamma*scaling_factor #compute scaled gamma
      Gamma_all[i] <- unique(Gamma[Gamma!=0]) #store scaled gamma
      
      squaresigma_y <- 1-PVE_zy
      
      epsilon_y <- rnorm(n2)
      epsilon_y <- epsilon_y/sd(epsilon_y)*sqrt(1-PVE_zy)
    
      y <- as.matrix(y_mean + epsilon_y)
      
      #run the PMR model using PMR_individual function
      yin=as.vector(x)
      zin=as.vector(y)
      x1in=zx
      x2in=zy
      result<-PMR_individual(yin, zin, x1in, x2in, method = "PMR_individual_Egger",  max_iterin = 1000, epsin = 1e-05,  Heritability_geneexpression_threshold = 1e-04)
      
      #get the estimate of the causal effect
      alpha_hat[i] <- result$causal_effect
      
      #get the estimate of the pleiotropy effect
      gamma_hat[i] <- result$pleiotropy_effect
      
      #get the pvalue for the causal test
      pvalue_alpha[i] <- result$causal_pvalue
      
      #get the pvalue for the pleiotropy effect
      pvalue_gamma[i] <- result$pleiotropy_pvalue
    }
    
    results[iter,] <- c(n_effects_zy,
                        PVE_zy,
                        mean(Gamma_all),
                        mean(gamma_hat),
                        mean(pvalue_alpha<0.05, na.rm=T))
    
  }
}

## The estimate of gene expression heritability explained by cis-SNPs is smaller than the threshold #### The pvalue of the causal effect test was assigned to be NA ##

results

   sparsity   pve        gamma    gamma_hat false_positive_rate
1        50 1e-04 0.0004361192 0.0004448700          0.05900000
2        50 1e-03 0.0013791301 0.0013727826          0.03800000
3        50 2e-03 0.0019503845 0.0019875726          0.04900000
4        50 1e-02 0.0043611924 0.0043928085          0.03700000
5         5 1e-04 0.0035777269 0.0003369432          0.04200000
6         5 1e-03 0.0111853841 0.0010993759          0.05900000
7         5 2e-03 0.0158495040 0.0015704796          0.07000000
8         5 1e-02 0.0353808651 0.0034478772          0.22100000
9         1 1e-04 0.0100000000 0.0002255487          0.05105105
10        1 1e-03 0.0316227766 0.0006625655          0.07500000
11        1 2e-03 0.0447213595 0.0009333676          0.11100000
12        1 1e-02 0.1000000000 0.0019830030          0.42400000

Conclusion

The false positive rate of PMR is well controlled under the dense scenario (sparsity=50), even when PVE of pleiotropy is large. This simulation matches the Eggar model underlying PMR. Holding PVE fixed, sparsity leads to inflated false positives. This inflation can be quite pronounced when the pleiotropy effects are very sparse (sparsity=1 or 5) with high PVE (PVE=0.01). Given that our cTWAS simulations are more similar to the very sparse scenarios, it is not surprising that PMR identifies many false positives in our simulations.

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] PMR_1.0

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.6         whisker_0.3-2      knitr_1.23        
 [4] magrittr_2.0.3     workflowr_1.6.2    R6_2.5.0          
 [7] rlang_1.0.2        fastmap_1.1.0      fansi_0.5.0       
[10] stringr_1.4.0      tools_3.6.1        xfun_0.8          
[13] utf8_1.2.1         cli_3.3.0          git2r_0.26.1      
[16] htmltools_0.5.2    ellipsis_0.3.2     rprojroot_2.0.2   
[19] yaml_2.2.0         digest_0.6.20      tibble_3.1.7      
[22] lifecycle_1.0.1    crayon_1.4.1       later_0.8.0       
[25] vctrs_0.4.1        fs_1.5.2           promises_1.0.1    
[28] glue_1.6.2         evaluate_0.14      rmarkdown_1.13    
[31] stringi_1.4.3      compiler_3.6.1     pillar_1.7.0      
[34] PDSCE_1.2.1        CompQuadForm_1.4.3 httpuv_1.5.1      
[37] pkgconfig_2.0.3